A general method to describe intersystem crossing dynamics in trajectory surface hopping

TL;DR

This paper reviews methods for modeling intersystem crossing, focusing on a new trajectory surface hopping approach called SHARC that includes spin-orbit and non-adiabatic couplings for large molecular systems.

Contribution

The paper introduces and details the SHARC method, a general trajectory surface hopping technique capable of incorporating arbitrary couplings, including spin-orbit interactions, for excited-state dynamics simulations.

Findings

SHARC can include any arbitrary coupling in surface hopping.

The method is applicable to large systems with many degrees of freedom.

Different electronic bases for spin-orbit coupling are analyzed.

Abstract

Intersystem crossing is a radiationless process that can take place in a molecule irradiated by UV-Vis light, thereby playing an important role in many environmental, biological and technological processes. This paper reviews different methods to describe intersystem crossing dynamics, paying attention to semiclassical trajectory theories, which are especially interesting because they can be applied to large systems with many degrees of freedom. In particular, a general trajectory surface hopping methodology recently developed by the authors, which is able to include non-adiabatic and spin-orbit couplings in excited-state dynamics simulations, is explained in detail. This method, termed SHARC, can in principle include any arbitrary coupling, what makes it generally applicable to photophysical and photochemical problems, also those including explicit laser fields. A step-by-step derivation of the main equations of motion employed in surface hopping based on the fewest-switches method of Tully, adapted for the inclusion of spin-orbit interactions, is provided. Special emphasis is put on describing the different possible choices of the electronic bases in which spin-orbit can be included in surface hopping, highlighting the advantages and inconsistencies of the different approaches.

Equations

$$m_A\frac{{\mathrm{d}}^2}{\mathrm{d}{t}^2}{\mathbf{R}}_A(t)=-{\nabla }_A{E}_{\beta }(\mathbf{R}(t)).$$

$${\mathbf{a}}_A(t)=$$

$${\mathbf{R}}_A(t+\Delta t)=$$

$${\mathbf{a}}_A(t+\Delta t)=$$

$${\mathbf{v}}_A(t+\Delta t)=$$

$$|\Psi ⟩=\sum _{\alpha }|{\psi }_{\alpha }⟩⟨{\psi }_{\alpha }|\Psi ⟩,$$

$$\text{i}\frac{\mathrm{d}}{\mathrm{d}t}|\Psi ⟩={\hat{H}}_{\text{el}}|\Psi ⟩$$

$$\text{i}\sum _{\alpha }⟨{\psi }_{\beta }|\frac{\mathrm{d}}{\mathrm{d}t}|{\psi }_{\alpha }⟩⟨{\psi }_{\alpha }|\Psi ⟩+⟨{\psi }_{\beta }|{\psi }_{\alpha }⟩\frac{\mathrm{d}}{\mathrm{d}t}⟨{\psi }_{\alpha }|\Psi ⟩=\sum _{\alpha }⟨{\psi }_{\beta }|{\hat{H}}_{\text{el}}|{\psi }_{\alpha }⟩⟨{\psi }_{\alpha }|\Psi ⟩.$$

$$\frac{\mathrm{d}}{\mathrm{d}t}⟨{\psi }_{\beta }|\Psi ⟩=-\sum _{\alpha }\left[\text{i}⟨{\psi }_{\beta }|{\hat{H}}_{\text{el}}|{\psi }_{\alpha }⟩+⟨{\psi }_{\beta }|\frac{\mathrm{d}}{\mathrm{d}t}|{\psi }_{\alpha }⟩\right]⟨{\psi }_{\alpha }|\Psi ⟩.$$

$$\frac{\mathrm{d}}{\mathrm{d}t}{\mathbf{c}}^{\text{rep}}=-\left[\text{i}{\mathbf{H}}^{\text{rep}}+{\mathbf{T}}^{\text{rep}}\right]{\mathbf{c}}^{\text{rep}}.$$

$${\hat{H}}_{\text{el}}^{\text{MCH}}=-\sum _i\left[\frac{1}{2}{\nabla }_i^2+\sum _A\frac{Z_A}{r_{Ai}}\right]+\sum _{i>j}\frac{1}{r_{ij}}+\sum _{A>B}\frac{Z_A{Z}_B}{r_{AB}}.$$

$${\hat{H}}_{\text{el}}^{\text{total}}={\hat{H}}_{\text{el}}^{\text{MCH}}+{\hat{H}}_{\text{el}}^{\text{add}},$$

$${\mathbf{H}}^{\text{diag}}={\mathbf{U}}^{†}{\mathbf{H}}^{\text{MCH}}\mathbf{U}.$$

$$H_{\beta \alpha }^{\text{diag}}={\delta }_{\beta \alpha }{E}_{\alpha }^{\text{diag}}.$$

$$\frac{\mathrm{d}}{\mathrm{d}t}{\mathbf{c}}^{\text{MCH}}=-\left[\text{i}{\mathbf{H}}^{\text{MCH}}+\mathbf{v}{\mathbf{K}}^{\text{MCH}}\right]{\mathbf{c}}^{\text{MCH}}.$$

$${\mathbf{c}}^{\text{diag}}={\mathbf{U}}^{†}{\mathbf{c}}^{\text{MCH}}.$$

$$\frac{\mathrm{d}}{\mathrm{d}t}{\mathbf{c}}^{\text{diag}}=-{\mathbf{U}}^{†}\left[\text{i}{\mathbf{H}}^{\text{MCH}}+\mathbf{v}{\mathbf{K}}^{\text{MCH}}\right]\mathbf{U}{\mathbf{c}}^{\text{diag}}-{\mathbf{U}}^{†}\frac{\mathrm{d}\mathbf{U}}{\mathrm{d}t}{\mathbf{c}}^{\text{diag}},$$

$${\mathbf{c}}^{\text{diag}}(t+\Delta t)=\exp \left[-\left(\text{i}{\mathbf{U}}^{†}{\mathbf{H}}^{\text{MCH}}\mathbf{U}+\mathbf{v}{\mathbf{U}}^{†}{\mathbf{K}}^{\text{MCH}}\mathbf{U}+{\mathbf{U}}^{†}\frac{\mathrm{d}\mathbf{U}}{\mathrm{d}t}\right)\Delta t\right]{\mathbf{c}}^{\text{diag}}(t),$$

$${\mathbf{c}}^{\text{MCH}}(t+\Delta t)=\exp \left[-\left(\text{i}{\mathbf{H}}^{\text{MCH}}+\mathbf{v}{\mathbf{K}}^{\text{MCH}}\right)\Delta t\right]{\mathbf{c}}^{\text{MCH}}(t),$$

$${\mathbf{c}}^{\text{diag}}(t+\Delta t)={\mathbf{U}}^{†}(t+\Delta t)\exp \left[-\left(\text{i}{\mathbf{H}}^{\text{MCH}}+\mathbf{v}{\mathbf{K}}^{\text{MCH}}\right)\Delta t\right]\mathbf{U}(t){\mathbf{c}}^{\text{diag}}(t).$$

$${\mathbf{c}}^{\text{diag}}(t+\Delta t)=\underset{{\mathbf{P}}^{\text{diag}}(t+\Delta t,t)}{\underbrace{{\mathbf{U}}^{†}(t+\Delta t){\mathbf{P}}^{\text{MCH}}(t+\Delta t,t)\mathbf{U}(t)}}{\mathbf{c}}^{\text{diag}}(t),$$

$${\mathbf{H}}^{\text{MCH}}(x)=\left(\begin{array}{cc}0.1x^2& \xi \\ \xi & 0.1(x-2{)}^2\end{array}\right).$$

$$h_{\beta \to \alpha }^{\text{diag}}=\left(1-\frac{{\left|c_{\beta }^{\text{diag}}(t+\Delta t)\right|}^2}{{\left|c_{\beta }^{\text{diag}}(t)\right|}^2}\right)\times \frac{\Re \left[c_{\alpha }^{\text{diag}}(t+\Delta t){\left(P_{\alpha \beta }^{\text{diag}}\right)}^{\ast }{\left(c_{\beta }^{\text{diag}}(t)\right)}^{\ast }\right]}{{\left|c_{\beta }^{\text{diag}}(t)\right|}^2-\Re \left[c_{\beta }^{\text{diag}}(t+\Delta t){\left(P_{\beta \beta }^{\text{diag}}\right)}^{\ast }{\left(c_{\beta }^{\text{diag}}(t)\right)}^{\ast }\right]}.$$

$$\sum _{i=1}^{\alpha -1}{h}_{\beta \to i}<r\le \sum _{i=1}^{\alpha -1}{h}_{\beta \to i}+h_{\beta \to \alpha }$$

$${\nabla }_{\mathbf{R}}⟨{\psi }_{\beta }^{\text{diag}}\left|{\hat{H}}_{\text{el}}^{\text{total}}\right|{\psi }_{\alpha }^{\text{diag}}⟩={\nabla }_{\mathbf{R}}{H}_{\beta \beta }^{\text{diag}}{\delta }_{\beta \alpha }.$$

$$⟨{\psi }_{\beta }^{\text{diag}}\left|{\nabla }_{\mathbf{R}}{\hat{H}}_{\text{el}}^{\text{total}}\right|{\psi }_{\alpha }^{\text{diag}}⟩={\nabla }_{\mathbf{R}}{H}_{\beta \beta }^{\text{diag}}{\delta }_{\beta \alpha }-(H_{\beta \beta }^{\text{diag}}-H_{\alpha \alpha }^{\text{diag}}){\mathbf{K}}_{\beta \alpha }^{\text{diag}},$$

$$⟨{\psi }_{\beta }^{\text{diag}}\left|{\nabla }_{\mathbf{R}}{\hat{H}}_{\text{el}}^{\text{total}}\right|{\psi }_{\alpha }^{\text{diag}}⟩=\sum _i\sum _j{U}_{i\beta }^{\ast }{U}_{j\alpha }⟨{\psi }_i^{\text{MCH}}|{\nabla }_{\mathbf{R}}({\hat{H}}_{\text{el}}^{\text{MCH}}+{\hat{H}}_{\text{el}}^{\text{add}})|{\psi }_j^{\text{MCH}}⟩.$$

$$⟨{\psi }_{\beta }^{\text{diag}}\left|{\nabla }_{\mathbf{R}}{\hat{H}}_{\text{el}}^{\text{total}}\right|{\psi }_{\alpha }^{\text{diag}}⟩=\sum _i\sum _j{U}_{i\beta }^{\ast }{U}_{j\alpha }⟨{\psi }_i^{\text{MCH}}\left|{\nabla }_{\mathbf{R}}{\hat{H}}_{\text{el}}^{\text{MCH}}\right|{\psi }_j^{\text{MCH}}⟩.$$

$${\nabla }_{\mathbf{R}}{H}_{\beta }^{\text{diag}}{\delta }_{\beta \alpha }-(H_{\beta \beta }^{\text{diag}}-H_{\alpha \alpha }^{\text{diag}}){\mathbf{K}}_{\beta \alpha }^{\text{diag}}=\sum _i\sum _j{U}_{i\beta }^{\ast }{U}_{j\alpha }\left[{\nabla }_{\mathbf{R}}{H}_i^{\text{MCH}}{\delta }_{ij}-(H_{ii}^{\text{MCH}}-H_{jj}^{\text{MCH}}){\mathbf{K}}_{ij}^{\text{MCH}}\right].$$

$$G_{ij}^{\text{rep}}={\nabla }_{\mathbf{R}}{H}_{ii}^{\text{rep}}{\delta }_{ij}-(H_{ii}^{\text{rep}}-H_{jj}^{\text{rep}}){\mathbf{K}}_{ij}^{\text{rep}},$$

$${\mathbf{G}}^{\text{diag}}={\mathbf{U}}^{†}{\mathbf{G}}^{\text{MCH}}\mathbf{U}.$$

$${\mathbf{P}}^{\text{MCH}}(t+\Delta t,t)=$$

Full text

Introduction

After a molecule is irradiated with light from the visible or UV parts of the spectrum, several photophysical processes can take place.1, 2, 3 These can be classified into radiative and radiationless processes and are best represented in a Jabłonski diagram (see Figure 1). In radiative processes, the molecule interacts with light and absorbs (absorption) or emits (stimulated and spontaneous emission) a photon. Emission can be further subdivided into fluorescence and phosphorescence, depending on whether the system emits from an electronic singlet or triplet state, respectively. Radiationless processes are usually divided into internal conversion (IC), which is the transfer of population between electronic states of the same spin multiplicity, and intersystem crossing (ISC), which involves a change of the electronic spin. Oftentimes, these processes occur in cascades after the incidental absorption of a photon. One typical chain of processes is absorption followed by IC to the lowest singlet state (as in Kasha’s rule4), followed by ISC to the triplet manifold, followed by IC to the lowest triplet state, from where the system may return to the electronic ground state by phosphorescence.

The non-radiative processes of IC and ISC are of fundamental importance in photochemistry and photobiology. IC is ubiquitous in the non-adiabatic excited-state dynamics of polyatomic molecules because polyatomic molecules possess many degrees of freedom and a large number of electronic excited states, allowing for regions of degeneracy or close degeneracy in their potential energy surfaces (PES). Intersections of two or more PES are usually refered to as conical intersections;5 these conical intersections are precisely the funnels which facilitate efficient non-adiabatic transitions from an upper to a lower PES, which is IC. IC is, for example, encountered in the process of visual perception, in particular during phototransduction. During this process, the retinal chromophore of the opsin protein absorbs a photon and subsequently undergoes ultrafast IC facilitated by conical intersections, leading to photoisomerization6, 7 which changes the conformation of the protein and induces a neural signal. Ultrafast IC is also responsible for avoiding harmful photochemical reactions of DNA8, 9 and proteins,10, 11, 12, 13, protecting organisms from UV radiation. IC is also present in the photochemistry of small molecules, for example O2 and O3,14 SO2,15 NO216 and other nitrogen oxides.17 The absorption of sunlight by these molecules and the subsequent photochemistry are central to atmospheric chemistry and the study of the environment and pollution.14

ISC is a key ingredient in molecular photophysics1, 3 and hence it plays an important role in many areas of research. For example, concepts like optical switching,18 reverse saturable absorption19 or light-induced magnetism20 are related to ISC in the field of molecular electronics. The interconversion between electric energy and light (in artificial photosynthesis21, 22 and OLEDs23), magnetic data storage20 and a number of biomimetic technologies24, 25 also rely on ISC. ISC has also biochemical relevance as it is a key step in many biochemical reactions (e.g., oxygen binding to carrier proteins 26), harmful photochemical reactions of DNA strands, and has found application in phototherapies (e.g., photodynamical therapy27). Furthermore, ISC is relevant to the study of combustion processes, where open-shell atoms (like oxygen) and radicals play a central role, in organic photochemistry, e.g., in $[2+2]$-cycloadditions,28 and in biomimetic catalysis.29

Due to its increasing relevance in different disciplines, the number of studies on ISC has proliferated in the last years. The study of photo-induced processes in transition-metal complexes,30 where ISC should be prominent due to the presence of the heavy metal atoms, has received special attention. Chergui31 has recently reviewed the photophysics of different transition metal complexes, paying particular attention to the timescale of ISC. Intriguingly, the reported ISC timescales range from extremely fast, e.g., 30 fs for the $[\text{Fe}(\text{bpy})_{3}]^{2+}$ complex,32 to about 800 ps in $\text{Pd}(\text{thpy})_{2}$.33 This large spread is difficult to explain simply by the type of central atom, since for example the heavier Pd is expected to induce more efficient ISC than Fe. Femtosecond time scales for ISC have been also observed in transition metal complexes based on ruthenium,34 rhenium,35 osmium,36 and other metals.31 More recently, also the role of ISC in organic molecules has captured the attention of researchers. Examples of molecules where ISC is important range from aldehydes37 and small aromatic compounds like benzene,38 naphthalene, anthracene and their carbonylic derivatives39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51 to nitrocompounds.43, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61 Furthermore, ISC has been reported for thio-substituted,62, 63, 64, 65, 66, 67, 68 aza-substituted,69 bromo-substituted70 and canonical nucleobases.71, 72, 73

From a theoretical point of view, ISC can be explained by the interaction of states of different multiplicity by spin-orbit coupling (SOC), which is a relativistic effect.74 It arises from the interaction of the magnetic field induced by the orbital momentum of an electron and the magnetic dipole related to the intrinsic electron spin. Since the orbital momentum depends on the nuclear charges, SOC is stronger in molecules possessing heavy atoms. For systems lacking heavy atoms, SOC has been historically regarded as negligible and hence ISC was believed to show much smaller reaction rates than IC.2 However, the actual ISC rate depends not only on the size of the SOCs, but also on the density of states, and the particular shape of the excited-state potential energy landscapes, most importantly the location and accessibility of surface crossing regions (e.g., between singlets and triplets). This is the reason why an a priori assessment of ISC rates is difficult and experiments have found that ISC can occur on very different time scales (from fs to ns), apparently not correlated with the magnitude of SOCs. In this context, theoretical methodologies are particularly useful to disclose the factors that determine ISC in a general framework.

For relatively slow ISC processes, the calculation of ISC rate constants using Landau-Zener theory,75 Fermi’s golden rule76, 77 or Marcus-Jortner theory76, 78, 79 is well established. For ultrafast ISC processes – where the system is very far from equilibrium – a dynamical treatment describing the motion on the PES and taking into account the kinetic processes that lead the system towards the relevant crossing region is mandatory. Encompassing the experimental discoveries involving ISC, the last decades have seen a surge of excited-state nuclear dynamics studies describing ISC. A number of studies have used quantum wavepacket dynamics, where inclusion of SOC is straightforward; however, these studies are usually restricted to few degrees of freedom. Examples include dihalogens in argon matrices,80, 81, 82 transition-metal complexes,83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95 and collision reactions.96 Further quantum dynamical studies investigated ISC in benzene,97 hydrogen fluoride,98 and sulphur dioxide.99 Full-dimensional dynamical studies can be performed using extensions of the trajectory surface hopping methodology. Different schemes have been employed to investigate ISC in collision reactions like $\text{O}+\text{C}_{2}\text{H}_{4}$,100, 101, 102, 103 $\text{O}+\text{N}_{2}$,104 $\text{Na}+\text{HCl}$,105 or $\text{S}+\text{H}_{2}$106 as well as in the dynamics of O2 on Al surfaces.107 Trajectory surface hopping was also used to treat ISC in molecules like sulphur dioxide,108 acrolein,109 acetone,110 pentanal,111 2-butene,112 6-thioguanine,113 cytosine,114, 115 uracil,116 as well as few transition-metal complexes.117, 118

The present contribution aims at describing and reviewing the available approaches described in the literature to treat ISC in trajectory surface hopping methods, paying special attention to the general surface hopping method recently developed by the authors,119, 120 which allows treating non-adiabatic IC and ISC on the same footing. This methodology is especially well suited to investigate ultrafast ISC (on timescales up to a few ps), because long timescales are usually too expensive for nuclear dynamics approaches. In this regard, our methodology can be considered complementary to other non-dynamical treatments of ISC,76, 77, 78, 79 which are not well applicable to ultrafast processes, but work well for slower reactions. Our methodology has been recently implemented in a completely new Sharc surface hopping dynamics code. The Sharc program suite, which contains the dynamics code itself and a number of setup, interface and analysis tools, is available online free of charge.120

The rest of the paper is then organized as follows. The next section introduces generally the surface hopping method and describes different ways in which SOC has been introduced to account for ISC. Then, section 3 describes step-by-step a way how to arrive at a generally applicable, computationally feasible and satisfactorily accurate methodology for the description of excited-state molecular dynamics including ISC. Finally, section 4 summarizes the presented methodology and gives an outlook on future developments in this field.

Dynamics Simulations of Intersystem Crossing

The basic equation for excited-state dynamics simulations is the time-dependent Schrödinger equation, where the Hamiltonian includes all non-adiabatic and potential couplings between the electronic states of interest. Simulating the motion of a quantum wavepacket according to this Hamiltonian by numerically integrating the time-dependent Schrödinger equation is usually termed quantum dynamics.121, 122 This approach includes all quantum-mechanical effects like tunneling, coherence and the splitting of the wavepacket and can deliver quantitative results given accurate PES.121 However, the necessity to know a priori the full multi-dimensional PES is a severe bottleneck of the method, rendering quantum dynamics in full dimensionality unfeasible for systems with many degrees of freedom. There exist a number of approximations that allow treating an increasing number of degrees of freedom, e.g. the “multi-configurational time-dependent Hartree” method (MCDTH),121 its variants “Gaussian-based MCTDH”,123 “multilayer MCTDH”124 and “variational multi-configurational Gaussians”,125 and the unrelated “full multiple spawning” method.126 Alternatively, the dimensionality problem can be tackled with trajectory surface hopping non-adiabatic dynamics.

The concept of surface hopping was originally devised by Tully127 and greatly improved later by his “fewest-switches criterion”.128 Given its wide acceptance and broad applicability, it has been reviewed extensively since then, see for example references129, 130, 131, 132. Surface-hopping approximates the motion of the nuclear wavepacket on the excited-state PES by an ensemble of a large number of independent, classical trajectories. The nuclear motion is determined by Newton’s equations, where the force acting on the nuclei is the gradient of one particular electronic state (the active state). Non-adiabatic population transfer between the electronic states can lead to instantaneous switches of the active state. This switching – the eponymous “hopping” – is conducted stochastically based on the rate of change of the electronic populations128 and is accompanied by an adjustment of the kinetic energy in order to preserve the total energy.

Surface-hopping methods can deliver a wealth of information about the excited-state dynamics of molecules. The ensemble ansatz allows obtaining branching ratios between different reaction channels. Time-dependent excited-state populations indicate the evolution of the system across different states, enabling to propose kinetic models. The geometries at which surface hops occurred can be used to infer surface crossing points, both conical intersections and singlet-triplet crossings (or crossing points between states of any multiplicity). By calculating observables along with the trajectories and taking the ensemble average, it is even possible to directly correlate the results with experimental observations.

Compared to other dynamics methodologies, surface hopping has a number of advantages,130 which justify its popularity: (i) The classical mechanics ansatz for the nuclear motion makes the method conceptually simple and easy to implement. (ii) The trajectories need only local information on the electronic states (i.e., energies, gradients, couplings). Thus, surface hopping allows for on-demand (“on-the-fly”) calculation of the electronic properties, which makes a treatment of all nuclear degrees of freedom possible without the need to construct the full multi-state PES beforehand. (iii) Each trajectory is independent from the rest of the ensemble. Hence, the trajectories can be trivially parallelized.

The fact that all degrees of freedom can be included in the simulation allows surface hopping to describe large systems, which is likely the reason of its widespread use. Another strength of on-the-fly methods is that they allow going beyond the so-called linear vibronic coupling133, 134 and linear vibronic spin-orbit coupling135, 77 approaches, where the relevant matrix elements are only approximated by some truncated power series around a reference geometry. In trajectory surface hopping methods, the availability of the full geometry dependence of the electronic properties makes it automatically possible to describe effects like the promotion of SOC by excitation of symmetry-reducing vibrational modes. Unfortunately, the classical approximation has the price that some quantum-mechanical effects are not treated properly. For example, nuclear vibration is not quantized, which means that nuclear interference effects are neglected and zero-point energy is formally not treated correctly. Additionally, because of the independent-trajectory ansatz, quantum coherence between the electronic states is not properly taken into account. The latter effect can be partially mitigated using so-called decoherence corrections; the interested reader is referred to the according literature.136, 137, 138, 139, 140, 141, 142, 143, 144, 145, 146, 147, 148, 149, 150, 151, 152, 153, 132 Tunneling effects can also not be described with surface hopping although there exist several approaches to alleviate this deficiency, see, e.g., references154, 155, 156, 157, 158.

The original surface hopping method was formulated to consider only non-adiabatic couplings128, 155 between electronic states of the same multiplicity and has been therefore intensively used to describe excited state dynamics involving IC. By incorporating SOC into the surface hopping procedure, dynamical studies can be extended to also describe ISC. In the following, we list the different concepts which – to our knowledge – have been described in the literature, ordered by how consistently ISC is described with regard to IC:

• •

On the simplest level, hopping between electronic states of different multiplicity is performed “manually”. This can be done by exploring first the deactivation mechanism within, e.g., the singlet manifold with regular surface hopping and then identifying singlet-triplet crossing regions, where ISC is assumed to occur. These regions then serve as a starting point for an independent simulation in the triplet manifold. This approach was used in the 70s by Warshel and Karplus,112 but is still in use nowadays for the description of exceptionally difficult systems.117, 111, 118

• •

On the next level, regular surface hopping can be combined with hopping probabilities between states of different multiplicities calculated only at crossing points based on Landau-Zener theory.159 Although in this way hoppings due to IC and hoppings due to ISC are treated inconsistently, this approach has been used to some extent in describing collision reactions.102, 104, 103, 105

• •

A further improvement is the proper inclusion of spin-orbit couplings in the equation of motion which governs the evolution of the electronic populations and hence the hopping probabilities. This allows in principle to have ISC in every time step of the simulation, not just at singlet-triplet crossings. In the simplest variant, the propagation is carried out using electronic states which are eigenstates of the total spin ($\hat{S}^{2}$) and spin projection ($\hat{S}_{z}$) operators. In this paper we term these states “MCH” states (referring to the “molecular Coulomb Hamiltonian”, see below for a detailed discussion), while in the literature they are also referred to as “spin-diabatic” states.160 Surface-hopping using these states has been adopted in several cases, either assuming a constant value for the SOC matrix elements100, 101 or with geometry-dependent SOC.107, 109, 161, 162

• •

Finally, the most advanced schemes to describe ISC employ a diagonalization of the electronic Hamiltonian including SOC, and basically conduct Tully’s surface hopping using the obtained eigenstates. Because the full electronic Hamiltonian is diagonal in the basis of these states, here we label it the “diagonal” basis (also called “spin-adiabatic” elsewhere160). The diagonal electronic basis offers a number of advantages, which will be discussed in detail below. This approach has been used already in 2004 to describe the $\text{S}+\text{H}_{2}$ reaction106 and recently saw a number of applications.110, 160, 113, 119, 114, 115, 108, 116 A general on-the-fly implementation of surface hopping in this fashion is available in the Sharc suite.119, 120

In passing we note that similar methodologies using the equivalents of the MCH163, 164 or diagonal165, 166, 119, 167, 168 electronic bases have also been developed in different contexts to describe other processes, like interactions of electromagnetic fields with molecular dipole moments.

The following section is devoted to explain exclusively the two last approaches to ISC surface hopping, since only these offer a consistent description of ISC versus IC.

Surface Hopping Methodology for Intersystem Crossing

In the following, we provide a step-by-step derivation of the main working equations necessary to perform trajectory surface hopping simulations focussing on the application to ISC problems. We start giving an overview over the workflow of a surface hopping simulation. Then, we develop the equation of motion for the electronic wavefunction, which will be our main concern here, since ISC arises from coupling of electronic states of different multiplicity. Afterwards, we discuss which electronic basis functions should be used in the propagation of the wavefunction, and the implications of this choice, in particular with respect to the integration of the equations of motion and the gradients used for the nuclear dynamics.

Surface Hopping

Surface-hopping is a special type of molecular dynamics (MD). In MD, the nuclear motion – which should be ideally described by quantum mechanics – is approximated with classical mechanics. The nuclear equation of motion is in this case given by Newton’s second law, which for the nuclear position $\mathbf{R}_{A}$ of atom $A$ reads as:

[TABLE]

Here, $m_{A}$ is the mass of atom $A$ and $-\nabla_{A}E_{\beta}(\mathbf{R}(t))$ is the gradient of the energy of the electronic state $\beta$ (we shall call this state the active state) with respect to the position of atom $A$; the energy depends on the geometry $\mathbf{R}$ of all atoms. In the standard formulation of MD, at each time step only the gradient $-\nabla_{A}E_{\beta}$ needs to be evaluated, e.g., by electronic structure methods, and performing the MD simulation consists of solely integrating equation (1). Integration from time $t$ to time $t+\Delta t$ can be effectively done with the velocity-Verlet algorithm,169 which can be expressed by the following equations:

[TABLE]

where $\mathbf{a}_{A}$ is the acceleration of atom $A$ and $\mathbf{v}_{A}$ is the velocity of atom $A$.

In order to incorporate non-adiabatic effects by means of the surface hopping scheme, a prescription of how to choose the correct active state $\beta$ at each time step is required. This choice is based on the evolution of an electronic wavefunction along the nuclear trajectory. Accordingly, for every time step the following steps need to be carried out:

For the current nuclear position $\mathbf{R}(t)$, the necessary electronic properties (e.g., electronic Hamiltonian and non-adiabatic couplings) are determined. 2. 2.

From these properties, the electronic wavefunction $\Psi(t)$ is propagated. 3. 3.

From the electronic wavefunction, the hopping probabilities $\mathbf{h}$ are calculated and the active state $\beta$ is determined stochastically. 4. 4.

The gradient of $E_{\beta}$ is used to calculate the new geometry $\mathbf{R}(t+\Delta t)$.

Based on this scheme, it can be seen that the evolution of $\Psi$ and $\mathbf{R}$ are coupled. In the following subsection, we shall focus on the evolution of $\Psi$ and the hopping probabilities $\mathbf{h}$ which can be derived from $\Psi$.

Electronic Equation-of-Motion

We expand the electronic wavefunction $|\Psi\rangle$ (we omit the explicit time-dependence for brevity) by using $\sum\limits_{\alpha}|\psi_{\alpha}\rangle\langle\psi_{\alpha}|=1$:

[TABLE]

where $\alpha$ runs over all electronic states in the model space. In order to arrive at a computationally tractable problem, one usually includes only a few states in the model space (e.g. the few lowest electronic states of each multiplicity for simulations of ISC dynamics), chosen in a way that the model space describes the total electronic wavefunction reasonably well for the processes under investigation.

In order to obtain the electronic equation of motion, we insert expression (6) into the electronic time-dependent Schrödinger equation,

[TABLE]

and project on $\langle\psi_{\beta}|$:

[TABLE]

Here, $\hat{H}_{\text{el}}$ is the electronic Hamiltonian operator, which we will define more precisely below. Note that in the following we refer to the electronic Hamiltonian simply by the “Hamiltonian”, as nuclei are treated classically and hence there is no nuclear Hamiltonian.

By using $\langle\psi_{\beta}|\psi_{\alpha}\rangle=\delta_{\beta\alpha}$, and after rearranging and dividing by i, one obtains:

[TABLE]

This equation is generally valid regardless the particular choice of the (complete) basis $|\psi_{\alpha}\rangle$, which we will call the “representation” in the following. The matrix representation of the Hamiltonian in a given basis “rep” will be denoted by bold $\mathbf{H}^{\text{rep}}$ in the remainder of the text, and the Hamiltonian matrix elements by $H_{\beta\alpha}^{\text{rep}}=\langle\psi_{\beta}^{\text{rep}}|\hat{H}_{\text{el}}|\psi_{\alpha}^{\text{rep}}\rangle$. For simplicity in the notation, the subscript “el” will by omitted in all electronic Hamiltonian matrices and matrix elements. We also define the wavefunction coefficient vector $\mathbf{c}^{\text{rep}}$ containing the elements $\langle\psi_{\alpha}^{\text{rep}}|\Psi\rangle$ and $\langle\psi_{\beta}^{\text{rep}}|\Psi\rangle$, and the temporal coupling matrix $\mathbf{T}^{\text{rep}}$ containing the elements $\langle\psi_{\beta}^{\text{rep}}|\frac{\mathrm{d}}{\mathrm{d}t}|\psi_{\alpha}^{\text{rep}}\rangle$. With these definitions, equation (9) can be written very compactly in the following matrix equation:

[TABLE]

Note that $\mathbf{T}^{\text{rep}}=\mathbf{v}\mathbf{K}^{\text{rep}}$, where $\mathbf{K}^{\text{rep}}$ is the non-adiabatic coupling matrix with elements $\langle\psi_{\beta}^{\text{rep}}|\nabla_{\mathbf{R}}|\psi_{\alpha}^{\text{rep}}\rangle$ and $\mathbf{v}$ is the nuclear velocity vector. The basis in which matrix representations of operators are expressed are given by a superscript of the matrix.

The reasons and implications of the choice of the representation of the electronic states are discussed in the following.

Representations

In a purely quantum-mechanical treatment of the nuclear motion the particular choice of the electronic basis does not influence the results, as long as the space spanned by the basis functions does not change. In contrast, due to its semiclassical nature, the results of surface hopping are not invariant to the choice of the electronic basis,131 as will be explained below.

For the purpose of this paper, it is useful to discuss three types of representations, in which $\mathbf{H}^{\text{rep}}$ and $\mathbf{K}^{\text{rep}}$ have different properties. We start by considering the molecular Coulomb Hamiltonian (MCH), which contains only the kinetic energy of the electrons and Coulombic interactions within the molecule, but neither external fields (e.g., electric fields) nor relativistic effects (e.g., SOC):

[TABLE]

Here, the indices $i$ and $j$ run over the electrons, $A$ and $B$ over the nuclei, $Z$ is the nuclear charge, $r$ is the distance between two particles, and atomic units have been used. The basis spanned by the (few lowest) eigenstates of the MCH Hamiltonian is denoted here as the MCH representation. In quantum chemistry, the MCH states are often called “adiabatic”, but here we shall refrain from using this term since it can lead to ambiguity once SOCs are introduced. Since the MCH is spin-independent and thus commutes with $\hat{S}^{2}$ and $\hat{S}_{z}$, the MCH states are also eigenstates of these operators; hence, they are usually labelled by their $S$ and $M_{S}$ values. Furthermore, because states of different $S$ or $M_{S}$ values are not coupled by the MCH, generally only a single multiplicity has to be considered in the dynamics (in most cases, singlet states due to their prevalence in closed-shell organic molecules). In the MCH representation, the matrix $\mathbf{H}^{\text{MCH,MCH}}$ (with matrix elements $H_{\beta\alpha}=\langle\psi_{\beta}^{\text{MCH}}|\hat{H}_{\text{el}}^{\text{MCH}}|\psi_{\alpha}^{\text{MCH}}\rangle$) is diagonal, but the elements of the non-adiabatic coupling matrix $\mathbf{K}^{\text{MCH}}$ are not zero (since the basis functions change with nuclear position) and may locally become very large. Since nearly all electronic structure codes have been implemented to calculate MCH states, this representation is the natural one to perform on-the-fly dynamics simulations.

Let us now consider additional terms in the Hamiltonian, in particular terms which are not included in the Hamiltonian employed by standard quantum chemistry software. These additional terms together with the MCH form the total Hamiltonian,

[TABLE]

where $\hat{H}_{\text{el}}^{\text{add}}$ represents the additional terms that are responsible for the process under investigation. These could be, for example, the interaction between molecular dipole moments and an external electric field $-\bm{\mu}\bm{\epsilon}^{\text{ext}}$, which allows describing absorption and stimulated emission of electromagnetic radiation. For the purpose of this contribution – the description of ISC – $\hat{H}_{\text{el}}^{\text{add}}$ is equal to the spin-orbit operator. Details on the calculation of SOC are beyond of the scope of this work and can be found elsewhere.170, 171, 172, 74 Upon inclusion of additional terms, the matrix representation of the total Hamiltonian in the MCH representation $\mathbf{H}^{\text{total,MCH}}$ (now with matrix elements $H_{\beta\alpha}=\langle\psi_{\beta}^{\text{MCH}}|\hat{H}_{\text{el}}^{\text{total}}|\psi_{\alpha}^{\text{MCH}}\rangle$) is not diagonal anymore, since the eigenstates of the MCH and the eigenstates of the total Hamiltonian are, in general, not identical. Importantly, the off-diagonal couplings in $\mathbf{H}^{\text{total,MCH}}$ are usually delocalized over the PES – a fact not optimal for surface hopping, since delocalized couplings can lead to a non-zero transition probability also far away from crossing regions. Such a situation necessitates a much larger number of trajectories to sample the process correctly because surface hops may occur in a much larger phase space volume.

Besides introducing off-diagonal elements in the Hamiltonian, the spin-orbit operator lifts the degeneracy between states with the same spatial wavefunction and same $S$ but different $M_{S}$ (the components of a multiplet). As a consequence, the investigation of ISC requires to take into account all components of the relevant spin multiplets explicitly and each multiplet of spin $S$ will contribute $2S+1$ states to the model space in the simulation. Importantly, in surface hopping simulations, the sum of the transition probabilities into all multiplet components should be independent of rotation of the molecule in the laboratory frame. As shown by Granucci et al.,160 this requirement is not fulfilled by surface hopping in the MCH representation (including SOCs as the non-diagonal elements of $\mathbf{H}^{\text{total,MCH}}$), which is a problem of the MCH representation in the presence of SOCs. In other words, if SOCs are included, the MCH basis should be regarded as sub-optimal for surface hopping.

A solution to this problem is to choose a different basis – one in which the total Hamiltonian matrix is diagonal. This can be achieved by a unitary transformation from the MCH basis to what we call the diagonal basis:

[TABLE]

Henceforth, we exclusively discuss matrix representations of the total Hamiltonian, so that $\mathbf{H}^{\text{rep}}$ always has matrix elements $H_{\beta\alpha}=\langle\psi_{\beta}^{\text{rep}}|\hat{H}_{\text{el}}^{\text{total}}|\psi_{\alpha}^{\text{rep}}\rangle$.

By definition, the matrix representation of the total Hamiltonian in the diagonal representation $\mathbf{H}^{\text{diag}}$ is diagonal.

[TABLE]

All couplings between the diagonal states are described by the non-adiabatic couplings $\mathbf{K}^{\text{diag}}$. Since the non-adiabatic couplings only become large when two PES are close to each other, the requirement that all couplings are localized is fulfilled. Additionally, the diagonalization solves the problem of rotational invariance of the multiplet components, since the eigenstates of the Hamiltonian do not depend on the molecular orientation: any rotation of the spin quantization axis could be described by a unitary transformation of the Hamiltonian $\mathbf{X}^{\dagger}\mathbf{H}\mathbf{X}$, which is similar to the original Hamiltonian and hence has the same eigenvalues. Since the diagonal representation (also refered to as ”fully adiabatic” in reference119) solves the problems of coupling localization and rotational invariance, and also best describes the energetics of the system, it can be regarded as the optimal representation to perform surface hopping for ISC.

For completeness, we also discuss the so-called “diabatic” representation.5 This representation is defined by having time- or geometry-independent basis functions (also called “crude adiabatic” basis173), which means that a diabatic basis has $\mathbf{T}^{\text{diab}}=0$ and likewise $\mathbf{K}^{\text{diab}}=0$. The coupling between diabatic states is described by off-diagonal terms in the Hamiltonian. This basis has the advantage that no diverging non-adiabatic couplings need to be considered, and the couplings in the Hamiltonian are easier to treat numerically. Another merit of using this electronic basis is that experimental observables are not strongly geometry-dependent and therefore, dynamical results in a diabatic basis facilitate discussion with respect to experiments. The disadvantage of this basis is that, in general, it is not clear whether strictly diabatic states for polyatomic molecules exist174, 175 and only a “quasi-diabatic” basis can be constructed, where most non-adiabatic couplings are small. An additional problem of the usage of a diabatic basis is the need to first diabatize the adiabatic states obtained from quantum chemistry, a procedure which is usually not possible with on-the-fly techniques. Hence, diabatic basis functions are in general not applicable to surface hopping simulations, at least not for full-dimensional treatments of larger molecules using the on-the-fly approach.

The three electronic representations discussed before are illustrated in Figure 2, which show examples for PES of a model system. On the left side, a two-state system is shown. In panel (a), two diabatic states, say ${}^{1}\pi\pi^{*}$ and ${}^{1}n\pi^{*}$, show a crossing in the center of the plot. For the purpose of this example, we assume an arbitrary total energy (indicated by a dashed line), which is not sufficient to reach the crossing and therefore access the right potential well. Hence, in a surface hopping simulation using these underlying potentials, the system would stay exclusively in the left minimum. In panel (b), the same states are plotted in the MCH representation. Now the electronic wavefunction character is not preserved and hence the states are called $S_{1}$ and $S_{2}$ (by energetic ordering). These states show a crossing strongly avoided, which lowers the energy barrier and makes accessing the right minimum classically allowed. In this representation, a classical dynamics simulation (e.g., with surface hopping), would visit both minima. The comparison of panels (a) and (b) clearly exemplifies the fact that the electronic representation has a direct influence on classical trajectories, which is one of the reasons why the results of surface hopping are indeed not invariant to the representation. Panel (c) shows that in the absence of additional couplings the MCH representation is identical with the diagonal representation. The distinction between MCH and diagonal representation comes into play only if $\hat{H}_{\text{el}}^{\text{add}}\neq 0$.

On the right side of figure 2, we show an exemplary three-state system (two singlets and one triplet), which are coupled by SOC. In panel (d), the diabatic PES of the ${}^{1}\pi\pi^{*}$ and ${}^{1}n\pi^{*}$ are unchanged compared to (a), since the SOCs only appear as off-diagonal elements of the Hamiltonian. The same holds true for panel (e), which is identical to (b) except for the addition of the triplet state. Only when transforming from the MCH to the diagonal representation (panel (f)) the PES are affected by SOC. Importantly, SOC lifts the degeneracy of the triplet components, giving rise to five states in total. A trajectory will explicitly depend on which triplet component it is moving in, again showing that the choice of the representation directly influences the result of classical trajectories, and that the diagonal representation performs best in this regard. Together with the fact that the diagonal representation optimally localizes all couplings in $\mathbf{H}^{\text{rep}}$ and $\mathbf{K}^{\text{rep}}$, we conclude that this representation is indeed optimal for surface hopping.

One should note, however, that the validity of other representations depends crucially on the strength of the coupling between the electronic states. If the SOCs are large, a large effect on the form of the PES can be expected and the diagonal representation is mandatory. For very weak SOC, the PES are only weakly deformed and other representations should yield approximately the same results, which justifies the usage of, e.g., the MCH representation in some forms of surface hopping including additional couplings.163, 164, 107, 109

Surface hopping in the diagonal representation

Currently there is no quantum chemistry program able to calculate all necessary properties needed for surface hopping in the diagonal basis. While gradients of diagonal states including SOC are already described for semi-empirical wavefunctions,176 no such gradients are available up to now for ab initio electronic structure methods and neither are the corresponding non-adiabatic couplings. Furthermore, quantum chemistry in the diagonal basis including SOCs (see references 171, 172, 177 and references therein) is in many cases significantly more expensive than quantum chemistry in the MCH basis. The higher cost is related to the fact that the spin-orbit Hamiltonian does not commute with $\hat{S}^{2}$ or $\hat{S}_{z}$ and hence it is not block-diagonal. Additionally, more demanding basis sets are usually required171, 172 for quantum chemical calculations in the diagonal basis than for MCH ones.

As a consequence, quantum chemistry in the MCH basis is nowadays the most advanced one. A number of quantum chemistry software suites have efficient implementations of analytical gradients, non-adiabatic couplings and SOC matrix elements,178, 179, 180 which can be employed for the on-the-fly calculations of surface hopping. Hence, a pragmatical approach for trajectory surface hopping methods is to combine quantum chemical data in the MCH representation with surface hopping in the diagonal representation. This is the essence of the ab initio molecular dynamics methodology Sharc,119 recently developed in our group. The basic idea is to construct the matrix representation of the total Hamiltonian in a small set of relevant MCH states (the model space) and diagonalize the Hamiltonian only in this small basis. The SOCs between the model space states and all remaining states are neglected, making this a form of quasi-degenerate perturbation theory (“QDPT”, see, e.g., reference 177). The obtained eigenstates are approximations to the true eigenstates of $\hat{H}_{\text{el}}^{\text{total}}$, but for small SOC this is a good approximation.

In the following, we describe how to perform surface hopping in a diagonal basis obtained from QDPT, as it is implemented in the Sharc code.120

Equation of motion

To make surface hopping in the diagonal representation practical, one needs a suitable equation of motion for the diagonal coefficients $\mathbf{c}^{\text{diag}}$ and a corresponding expression for the hopping probabilities $\mathbf{h}^{\text{diag}}$ based on the MCH quantities $\mathbf{H}^{\text{MCH}}$ and $\mathbf{K}^{\text{MCH}}$. Additionally, one needs to construct the diagonal gradient $-\nabla E_{\beta}^{\text{diag}}$ from the MCH gradients.

We start with the equation of motion (10) in the MCH basis:

[TABLE]

The diagonal wavefunction coefficients $\mathbf{c}^{\text{diag}}$ can be easily obtained from $\mathbf{c}^{\text{MCH}}$, using the transformation matrix $\mathbf{U}$:

[TABLE]

Inserting the latter into equation (15) and pre-multiplying with $\mathbf{U}^{\dagger}$, we obtain:

[TABLE]

Note how the last term on the right-hand side arises since the transformation matrix $\mathbf{U}$ implicitly depends on time via the time dependence of the Hamiltonian. Given a coefficient vector $\mathbf{c}^{\text{diag}}(t)$ and the matrices $\mathbf{H}^{\text{MCH}}$ and $\mathbf{K}^{\text{MCH}}$, equation (17) can be integrated using standard Runge-Kutta-type numerical methods, or using the exponential operator method for sufficiently short time steps:

[TABLE]

where we omitted the explicit time-dependence of $\mathbf{U}$, $\mathbf{H}^{\text{MCH}}$ and $\mathbf{K}^{\text{MCH}}$ for brevity. However, in this scheme the term $\mathbf{U}^{\dagger}\frac{\mathrm{d}\mathbf{U}}{\mathrm{d}t}$ leads to serious numerical difficulties. First, the unitary matrix $\mathbf{U}$ is not uniquely defined by equation (13) – each eigenvector could be multiplied by a complex phase factor $\text{e}^{\text{i}\phi}$ and still be an eigenvector to the same eigenvalue. This non-uniqueness makes it impossible to calculate $\mathbf{U}^{\dagger}\frac{\mathrm{d}\mathbf{U}}{\mathrm{d}t}$ from finite differences, because $\mathbf{U}$ is not continuous. Second, $\mathbf{U}^{\dagger}\frac{\mathrm{d}\mathbf{U}}{\mathrm{d}t}$ may become large if two MCH states cross and the norm of their mutual coupling tends to zero, which can generally occur in all systems. Using shorter time steps $\Delta t$ alleviates the problem, but the numerical precision of usual floating point calculations limits the size of the time step. Furthermore, if the mutual coupling is zero, $\mathbf{U}^{\dagger}\frac{\mathrm{d}\mathbf{U}}{\mathrm{d}t}$ diverges.

These problems can be easily circumvented by integrating (for small $\Delta t$) the equation of motion in the MCH basis (15),

[TABLE]

and only then inserting equation (16):

[TABLE]

In this way, one obtains a propagation equation which is fully equivalent to (18), but where $\mathbf{U}^{\dagger}\frac{\mathrm{d}\mathbf{U}}{\mathrm{d}t}$ is not explicitly involved. Note that this equation describes a propagation which proceeds in three distinct steps:

transforming $\mathbf{c}^{\text{diag}}(t)$ to $\mathbf{c}^{\text{MCH}}(t)$, 2. 2.

propagating $\mathbf{c}^{\text{MCH}}(t)$ to $\mathbf{c}^{\text{MCH}}(t+\Delta t)$, and 3. 3.

transforming $\mathbf{c}^{\text{MCH}}(t+\Delta t)$ to $\mathbf{c}^{\text{diag}}(t+\Delta t)$.

We call this approach the three-step propagator, and write it in short as

[TABLE]

where $\mathbf{P}^{\text{MCH}}(t+\Delta t,t)$ is the propagator matrix in the MCH basis, which acts on the MCH coefficients at time $t$ to propagate it to time $t+\Delta t$. We note here that $\mathbf{P}^{\text{MCH}}(t+\Delta t,t)$ could also be calculated differently from equation (20), for example using the local diabatization procedure of Granucci et al.181, 182 The condition that the time step $\Delta t$ is short enough can be fulfilled by choosing a shorter time step for the electronic propagation than for the nuclear propagation and linearly interpolating the electronic properties (see Appendix A for the corresponding equations).

The three-step-propagator methodology has a number of advantageous properties. First, it accomplishes our main goal, i.e., formulating the propagation of the diagonal coefficients based on MCH quantities. Second, it is unitary (the matrix exponential of an anti-Hermitian matrix is always unitary), and thus the wavefunction norm is automatically conserved during the propagation. And third, the numerical computation of $\mathbf{U}^{\dagger}\frac{\mathrm{d}\mathbf{U}}{\mathrm{d}t}$ is avoided, providing a much more stable propagation with respect to the time step and the involved couplings in the Hamiltonian.

The following example shows the numerical superiority of the three-step propagator over the more naïve single-step method in equation (18). We consider a simple one-dimensional system consisting of two harmonic oscillators coupled by a constant off-diagonal term $\xi$ (atomic units):

[TABLE]

The initial position of the system is set at $x=10$ in the upper diagonal state and the reduced mass is $m=0.2$. We simulate a single pass of a trajectory through the crossing region and record the total change in the electronic populations depending on the strength of $\xi$. We performed the simulations with equation (20) and a time step of $10^{-2}$ fs. For comparison we also used equation (18), with different time steps between $10^{-4}$ and $10^{-8}$ fs. In order to show that controlling the complex phases of the columns of $\mathbf{U}$ is imperative when using equation (18), results with and without phase control (using the algorithm outlined in Appendix B) are shown. Figure 3 shows the population of the initially empty state ($|c^{\text{MCH}}_{1}(0)|^{2}=0$) after passing the crossing. According to Landau-Zener theory, for sufficiently small couplings, the population after the crossing $|c^{\text{MCH}}_{1}(t>t_{\text{cross}})|^{2}$ should be proportional to $\xi^{2}$.

The three-step propagator shows the expected proportionality of $|c^{\text{MCH}}_{1}(t>t_{\text{cross}})|^{2}$ and $\xi^{2}$, even though a comparatively long time step was employed. Even if $\xi$ is set equal to zero the three-step propagator delivers the correct result of zero population transfer. The one-step propagator data shows that equation (18) only performs well when the coupling $\xi$ is large, while for smaller couplings, population transfer is systematically overestimated. The use of comparatively smaller time steps helps in the case of small couplings, but even the shortest time step of $10^{-8}$ fs fails for $\xi<10^{-3}$ cm*-1* when using the one-step propagator. When $\xi$ tends to zero, equation (18) predicts that the population $|c^{\text{MCH}}_{1}(t>t_{\text{cross}})|^{2}$ tends to one, which is of course unphysical. Not controlling the phase of $\mathbf{U}$ leads to a completely erratic behaviour and a much too large population transfer, which is basically independent from $\xi$.

The message from this figure is that the one-step propagator is very dangerous to use in the presence of very small couplings. This is a very serious problem in ISC simulations because, even though a limited number of SOC matrix elements may be large, the problematic small SOC matrix elements are generally present in all systems. Especially in small molecules, the selection rules of SOC render some matrix elements equal to zero, and the El-Sayed rule183 predicts that SOCs between states with very similar wavefunction character are small. Therefore, the use of the three-step propagator method as outlined above is recommended to avoid the problems in the wavefunction propagation in such situations.

Hopping probabilities

The use of the matrix exponential for the propagation of the coefficients is not compatible with the original expression for the hopping probabilities given by Tully.128 In the diagonal representation, this expression contains $\frac{\mathrm{d}}{\mathrm{d}t}\mathbf{c}^{\text{diag}}$, which according to equation 17 explicitly depends on $\mathbf{U}^{\dagger}\frac{\mathrm{d}\mathbf{U}}{\mathrm{d}t}$, giving rise to the numerical problems mentioned above. Consequently, it is necessary to obtain a different expression for the hopping probabilities $h_{\beta\rightarrow\alpha}$ – yet based on the fewest switches criterion. Here we use a variant of the equation derived in reference 181, where the hopping probability $h_{\beta\rightarrow\alpha}$ is just a function of the old and new coefficients and the propagator matrix:

[TABLE]

Negative hopping probabilities are set to zero, as is $h_{\beta\rightarrow\beta}$. The equation can be derived similarly as in reference181, but using $\mathbf{P}^{\dagger}\mathbf{c}(t+\Delta t)=\mathbf{c}(t)$ instead of $\mathbf{c}(t+\Delta t)=\mathbf{P}\mathbf{c}(t)$. One advantage of equation (23) over the one in reference181 is that the denominator in equation (23) does not tend to zero when the population transfer becomes large, giving more accurate results. Moreover, equation (23) has the nice property that the hopping probabilities are not affected by the non-uniqueness of the transformation matrix $\mathbf{U}$. Hence, no phase adjustment between $\mathbf{U}(t)$ and $\mathbf{U}(t+\Delta t)$ is necessary at all in the Sharc methodology.

After the propagation for a particular time step has been performed, equation (23) can be directly used to calculate hopping probabilities. To choose the active state for the next time step, a random number $r$ from the interval $0\leq r\leq 1$ is selected. A hop to state $\alpha$ is performed, if

[TABLE]

i.e., if the random number $r$ lies in an interval with a width proportional to the hopping probability.

Gradients

Equations (13), (20) and (23) define all that is needed to perform the electronic part of the surface hopping algorithm, namely to determine the active state for the next time step. The remaining necessary element of an algorithm for surface hopping in the diagonal representation is the construction of the gradient of the diagonal states.

We start with the derivative of a matrix element of the total Hamiltonian in the diagonal basis:

[TABLE]

By applying the product rule to the left-hand side and using the fact that the $|\psi^{\text{diag}}\rangle$ are the eigenfunctions of $\hat{H}_{\text{el}}^{\text{total}}$, equation (25) can be rearranged to129

[TABLE]

where $H_{\beta\beta}^{\text{diag}}$ and $H_{\alpha\alpha}^{\text{diag}}$ are eigenvalues of the total Hamiltonian and $\mathbf{K}_{\beta\alpha}^{\text{diag}}$ is a matrix element of the non-adiabatic coupling matrix $\mathbf{K}^{\text{diag}}$ (note that the elements of $\mathbf{K}^{\text{diag}}$ are itself vectors in the space of the nuclear distortions). Importantly, this equation contains $\nabla_{\mathbf{R}}H_{\beta\beta}^{\text{diag}}$, which is the diagonal gradient necessary for the dynamics.

The left-hand side of equation (26) can be written in terms of the MCH states and the transformation matrix $\mathbf{U}$ as

[TABLE]

The critical term in this equation is the derivative of the additional Hamiltonian terms with respect to the nuclear displacements $\nabla_{\mathbf{R}}\hat{H}_{\text{el}}^{\text{add}}$. If – as in the present case – $\hat{H}_{\text{el}}^{\text{add}}$ is a spin-orbit operator, then this derivative is (up to our knowledge) not available from standard ab initio software. However, due to the short-range nature of the spin-orbit interaction,184 the derivative $\nabla_{\mathbf{R}}\hat{H}_{\text{el}}^{\text{add}}$ is expected to be small, and is hence ignored. The remainder of equation (27) can thus be written as

[TABLE]

As the $\psi_{i}^{\text{MCH}}$ are the eigenfunctions of $\hat{H}_{\text{el}}^{\text{MCH}}$, we can use an analogue form of equation (26) to arrive at

[TABLE]

If we define the general gradient matrix in an arbitrary representation $\mathbf{G}^{\text{rep}}$ with matrix elements

[TABLE]

then we can write the above equation compactly as

[TABLE]

Equations (30) and (31) give a rule how to obtain the diagonal gradients and non-adiabatic couplings from the MCH gradients and non-adiabatic couplings. From an algorithmic point of view, the diagonal vector properties can be obtained in the following steps:

Obtain the energies $H_{ii}^{\text{MCH}}$, gradients $\nabla_{\mathbf{R}}H_{ii}^{\text{MCH}}$ and non-adiabatic couplings $\mathbf{K}_{ij}^{\text{MCH}}$ from quantum chemistry. 2. 2.

Construct $\mathbf{G}^{\text{MCH}}$ from these quantities using equation (30). 3. 3.

Transform $\mathbf{G}^{\text{MCH}}$ into the diagonal basis using equation (31). 4. 4.

Obtain the diagonal gradients $\nabla_{\mathbf{R}}H_{\beta\beta}^{\text{diag}}=\mathbf{G}^{\text{diag}}_{\beta\beta}$. 5. 5.

If needed, obtain the non-adiabatic coupling vectors in the diagonal basis $\mathbf{K}_{\beta\alpha}^{\text{diag}}=(H_{\beta\beta}^{\text{diag}}-H_{\alpha\alpha}^{\text{diag}})^{-1}\mathbf{G}^{\text{diag}}_{\beta\alpha}$.

As one can see from the last paragraphs, the gradients in the diagonal basis are linear combinations of the MCH gradients and the non-adiabatic coupling vectors. Therefore, even if only a single diagonal gradient is needed, several MCH gradients and coupling vectors need to be obtained from quantum chemistry. Compared to dynamics in the MCH basis, where only a single gradient needs to be computed at each time step, for dynamics in the diagonal basis, in principle, all gradients are required, adding computational cost to the simulation.

As a consequence of the above-said, the surface hopping procedure needs the non-adiabatic coupling vectors for two distinct steps: the wavefunction propagation and the gradient transformation. Unfortunately, the calculation of non-adiabatic couplings is presently only available in some quantum chemistry codes. If non-adiabatic couplings are not available, it is still possible to calculate the propagator matrix $\mathbf{P}$ using other schemes, see, e.g., references 181, 155. However, for the gradient transformation no alternative schemes without non-adiabatic couplings are available, and, as an approximation $\mathbf{K}_{ij}^{\text{MCH}}$ has to be neglected in equations (30) and (31). This introduces an error which is largest near weakly-avoided crossings (precisely where $\mathbf{K}_{ij}^{\text{MCH}}$ is large). Figure 4 shows the total energy of two exemplary trajectories,111The simulations were performed on seleno-acrolein at the CASSCF(6,5)/ano-rcc-vdzp level of theory including two singlet and two triplet electronic states. The initial geometry was the ground state minimum geometry and the initial state the $S_{1}$. The initial velocities of all atoms were zero, except for the terminal hydrogen atoms, which had velocities of 0.008 a.u. out of the molecular plane and in opposite directions, thus giving some initial momentum for torsion around the C=C double bond.177 one calculated with equations (30) and (31) and another where $\mathbf{K}_{ij}^{\text{MCH}}$ is neglected. As it can be seen in this extreme case, neglecting the contribution of the MCH non-adiabatic coupling vectors to the diagonal gradients leads to a violation of energy conservation already after a few fs. Hence, whenever possible it is recommended to include the non-adiabatic couplings in the gradient calculations or, if these are not included, one should check the total energy conservation carefully.

An algorithm for surface hopping in the diagonal representation

To summarize, given some initial geometry $\mathbf{R}(t)$, velocity $\mathbf{v}(t)$, acceleration $\mathbf{a}(t)$, coefficients $\mathbf{c}(t)$ and active state $\beta(t)$, the following steps should be performed to propagate the dynamics to time $t+\Delta t$:

Use equations (2) and (3) to find the new geometry $\mathbf{R}(t+\Delta t)$. 2. 2.

Obtain $\mathbf{H}^{\text{MCH}}$, $\mathbf{K}^{\text{MCH}}$ and the gradients $-\nabla H^{\text{MCH}}_{ii}$ from quantum chemistry software. 3. 3.

Find $\mathbf{U}$ from equation (13). 4. 4.

Use equations (30) and (31) to find $-\nabla H^{\text{diag}}_{\beta\beta}(t+\Delta t)$. 5. 5.

Use equations (4) and (5) to find the new velocity $\mathbf{v}(t+\Delta t)$. 6. 6.

Use equation (20) to find the propagator $\mathbf{P}(t+\Delta t,t)$ and the new coefficients $\mathbf{c}(t+\Delta t)$. 7. 7.

Use equation (23) to find the hopping probabilities $h_{\beta\rightarrow\alpha}$ and equation (24) to find the new active state $\beta(t+\Delta t)$. 8. 8.

If a hop occured, ensure total energy conservation by adjusting $\mathbf{v}(t+\Delta t)$ as in regular surface hopping (see, e.g., references 127, 185). 9. 9.

Apply some form of decoherence correction140, 136, 137 to the diagonal coefficients $\mathbf{c}(t+\Delta t)$.

A few comments regarding the different electronic state representations encountered in this workflow are appropriate at this point. First, it should be mentioned that using MCH quantum chemistry data to perform dynamics in the diagonal basis is actually an approximation based on a QDPT treatment and hence neglects the coupling to states outside the model space. Moreover, the fully consistent diagonal gradients could only be calculated if the derivatives of the SOCs were known. However, in view of the current state-of-the-art of quantum chemistry, dynamics directly in the diagonal basis is not possible because the necessary properties are not available. Second, while the diagonal representation is optimal for surface hopping dynamics, it is not optimal for the interpretation of the results, since the character and spin of the diagonal states is geometry-dependent. Therefore, it is usually advantageous to transform the results (e.g., populations) back into the MCH basis, where the states can be characterized by their $S$ and $M_{s}$ quantum numbers. This facilitates, for example, the distinction between IC and ISC processes, which are otherwise equivalent in the diagonal picture. Of course, this a posteriori transformation to MCH states is reasonable only for systems with SOC small enough so that the $LS$ coupling scheme is valid. Third, a transformation of the results from the MCH into a “spectroscopic” basis (an approximate diabatic basis where spectroscopic observables depend strongly on the state, but weakly on geometry, e.g., states of $\pi\pi^{*}$ and $n\pi^{*}$ character) is even more convenient than the MCH basis to facilitate the interpretation of the results and is recommended as long as an a posteriori diabatization scheme can be devised for the system.

Conclusion and Outlook

In this contribution it has been shown how spin-orbit coupling (SOC) can be treated on the same footing as non-adiabatic coupling in semi-classical surface hopping trajectory methods. To this aim, a general and accurate approach relying on the diagonalization of the electronic Hamiltonian including all potential couplings has been described. This methodology was coined “surface hopping including arbitrary couplings” (Sharc) in our original publication119 because it is not limited to SOC, but can treat all kinds of couplings in the same fashion. An on-the-fly implementation of this methodology is available in the Sharc molecular dynamics suite.120

While the Sharc methodology is generally applicable to very large systems (due to the classical nuclear dynamics approach), the size of the systems for which ISC can be studied is actually limited by the availability of suitable quantum chemistry methods. These methods need to provide accurate excitation energies, analytical gradients, SOC matrix elements, and non-adiabatic couplings (or equivalent properties like overlap matrices181, 182) for several excited states at once. If laser interactions are considered in the framework of Sharc (few applications are found in the literature so far119, 167, 168, 139), then the quantum chemistry also needs to provide transition dipole moments. Additionally, the quantum chemical methods should be able to give a balanced and accurate description of all states over the whole potential energy surface, which usually requires multi-configurational wavefunctions. These requirements are currently fulfilled for several implementations180, 179, 178 of the multi-configurational self-consistent field (MCSCF) method186, 187 and for multi-reference configuration interaction (MRCI) as implemented in Columbus.178 Using MCSCF or MRCI, already a number of ISC dynamics simulations have been reported using Sharc.114, 116, 115, 108, 188, 189

The progress of general ab initio molecular dynamics schemes like Sharc strongly depends on the development of electronic structure methods that can go beyond MCSCF, which is one of the most popular methods for on-the-fly excited-state dynamics. A obvious step further is complete-active-space perturbation theory (e.g., CASPT2190), which offers very good accuracy (oftentimes better than MRCI) and performance (usually much faster than MRCI). Few surface hopping studies based on CASPT2 have been reported already,191, 192, 193, 194 albeit no ISC dynamics simulation has been published until now, since no current implementation of CASPT2 can provide all necessary quantities. Other promising electronic structure methods are, e.g., the approximately size-consistent MRCI variant linear-response-theory multi-reference average-quadratic coupled-cluster (LRT-MRAQCC)177 or density matrix normalization group (DMRG) approaches.195 Furthermore, it will be possible to employ time-dependent density functional theory (TDDFT) in ISC dynamics studies with the help of recent developments,161, 162 although care has to be taken that the excited-state potential energy surfaces are described correctly in the case of strongly multi-configurational state character.

Semi-empirical electronic structure methods for excited states196 are very efficient in comparison to ab initio or TDDFT methods and can be applied as well to ISC dynamics simulations. There is already an implementation equivalent to Sharc,160 which has been applied to study ISC dynamics on thioguanine.113 Additionally, semi-empirical calculations of ISC dynamics already have been carried out using a conventional surface hopping approach,109 i.e., without the advantages provided by the Sharc methodology as outlined in this contribution. In any case, the availability of TDDFT and semi-empirical methods for surface hopping will open wide opportunities to describe larger molecules or achieve longer propagation times than it is currently feasible with multi-configurational wavefunction methods.

A further challenge in trajectory surface methods is to include an efficient and accurate description of solvents and complex molecular environments in the dynamics. Considerable effort in this direction has already been made for standard surface hopping using mixed quantum-mechanics/molecular-mechanics (QM/MM) methods, in which the main chromophore is treated quantum-chemically and the environment with force fields.197, 198, 199, 200

Finally, in order to relate the computational results to experimental data, not only the excited-state dynamics has to be simulated, but also the subsequent processes occuring during the experimental probing, i.e., it is necessary to actually simulate the experimental observables. A typical example is transient absorption spectroscopy,201 which can be simulated by including highly excited states and the corresponding transition dipole moments in the calculations. However, an accurate description of these highly excited states is far from trivial. Another typical probing technique is based on ionization of the studied molecule, yielding for example time-dependent photo-electron spectra.202 In order to simulate such ionizations, it is necessary to describe the wavefunctions of the neutral system, the ionized species as well as the outgoing electron; from these wavefunctions, the corresponding wavefunction overlaps as well as the associated transition dipole moments need to be calculated. This is quite a formidable task that has not yet been achieved in full, so that a detailed understanding of molecular ionization is still out of reach. Therefore, further developments are necessary to achieve the ultimate goal of simulating all the processes involved in actual experiments from the beginning to the end. Yet, given the increasing interest in dynamical processes including ISC, it is expected that the coming years will see an increasing number of surface hopping trajectory method simulations.

ACKNOWLEDGMENTS

This work is supported by the Austrian Science Fund (FWF) through the project P25827. The authors are specially thankful to the past and present members of the development team of Sharc – M. Richter, J. González-Vázquez, I. Sola, M. Oppel and M. Ruckenbauer – for their many contributions and discussions, and also to the many colleagues who through their questions and suggestions have contributed to improve Sharc from the first to the current implementation. We are also indebted to T. Müller, T. Hoffmann-Ostenhof, M. Persico and G. Granucci for inspiring discussions. The COST Actions CM1204 (XLIC) and CM1305 (ECOSTBio) are thanked for promoting different collaborative applications related to Sharc. Finally, the Vienna Scientific Cluster (VSC) is gratefully acknowledged for generous allocation of computer time.

Appendix A. Propagation using substeps

Equation (21) contains the propagation matrix $\mathbf{P}^{\text{MCH}}(t+\Delta t,t)$ in the MCH basis. For small $\Delta t$, this matrix can be calculated from equation (20). For too large $\Delta t$, it is possible to split the interval into $n$ subintervals $\Delta\tau=\frac{\Delta t}{n}$. Here, we give the corresponding equations for the calculation of $\mathbf{P}^{\text{MCH}}(t+\Delta t,t)$ using subintervals.

The propagator can be calculated by linearly interpolating $\mathbf{H}^{\text{MCH}}$ and $\mathbf{T}^{\text{MCH}}$, using

[TABLE]

where the matrix elements of $\mathbf{T}^{\text{MCH}}$ are

[TABLE]

(classical Tully scheme128) or

[TABLE]

(Hammes-Schiffer-Tully scheme155).

Using the local diabatization scheme,181, 182 the propagator is, instead, given by

[TABLE]

where $\mathbf{H}^{\text{MCH}}_{\text{tra}}$ is the matrix representation of $\hat{H}_{\text{el}}(t+\Delta t)$ in the basis of the states at time $t$. The overlap matrix $\mathbf{S}^{\text{MCH}}(t,t+\Delta t)$ has elements

[TABLE]

Appendix B. Continuous transformation-matrix phase

A Hermitian matrix $\mathbf{H}^{\text{MCH}}$ can always be diagonalized. Its eigenvectors form the rows of a unitary matrix $\mathbf{U}$, which can be used to transform between the original basis and the basis of the eigenfunctions of $\mathbf{H}$:

[TABLE]

However, the condition that $\mathbf{U}$ diagonalizes $\mathbf{H}^{\text{MCH}}$ is not sufficient to define $\mathbf{U}$ uniquely. Each normalized eigenvector can be multiplied by a complex number on the unit circle and still remains a normalized eigenvector. We define a second transformation matrix $\tilde{\mathbf{U}}=\mathbf{U}\bm{\Phi}$, where $\bm{\Phi}$ is a diagonal matrix with elements $\Phi_{\beta\alpha}=\delta_{\beta\alpha}\text{e}^{\text{i}\phi_{\alpha}}$ (i.e., an eigenvector of $\tilde{\mathbf{U}}$ is the same eigenvector of $\mathbf{U}$ times a complex phase). By this definition, $\bm{\Phi}$ is unitary ($\bm{\Phi}^{\dagger}\bm{\Phi}=\mathbf{1}$) and commutes with $\mathbf{H}^{\text{diag}}$, since two diagonal matrices always commute. With this properties it can be shown by

[TABLE]

that also $\tilde{\mathbf{U}}$ diagonalizes $\mathbf{H}^{\text{MCH}}$ and produces the same eigenvalues.

In most libraries for numerical diagonalization of Hermitian matrices, the eigenvector phases given in $\bm{\Phi}$ are either uncontrolled (so that the phases depend on the details of the implementation) or are fixed by some constraints (e.g., the first non-zero element of the eigenvector is made real-valued). However, by these conditions the eigenvector phases are not continuous in time – a small change in geometry from $\mathbf{R}(t)$ to $\mathbf{R}(t+\Delta t)$ may lead to a large jump in some eigenvector phases. Since the transformation matrix with discontinuous phases $\mathbf{U}(t)\bm{\Phi}(t)$ is itself discontinuous, it is not differentiable. Thus, it is not possible to calculate the antihermitian term $\mathbf{U}^{\dagger}\frac{\mathrm{d}}{\mathrm{d}t}\mathbf{U}$ from equation (17), neither analytically nor by finite differences. Attempting to calculate finite differences anyways will lead to a result that is not antihermitian and generally overestimates the size of the matrix elements of the differential. Thus, attempting to use equation (17) for the propagation of the wavefunction and relying on the finite-difference calculation of $\mathbf{U}^{\dagger}\frac{\mathrm{d}}{\mathrm{d}t}\mathbf{U}$ without making the eigenvector phases (approximately) continuous will lead both to overestimated population transfer – possibly at any time – and to violation of the conservation of the wavefunction norm.

In order to facilitate the calculation of $\mathbf{U}^{\dagger}\frac{\mathrm{d}}{\mathrm{d}t}\mathbf{U}$, it is necessary to control the phase of the transformation matrix $\mathbf{U}$. Specifically, this means that the phases of $\mathbf{U}(t+\Delta t)$ have to be aligned to the phases of $\mathbf{U}(t)$. Hence, given $\mathbf{U}(t)$ and $\tilde{\mathbf{U}}(t+\Delta t)$, we need to find $\bm{\Phi}$, which allows to obtain the phase-aligned transformation matrix $\mathbf{U}(t+\Delta t)$:

[TABLE]

Choosing $\bm{\Phi}^{\dagger}=\tilde{\mathbf{U}}^{\dagger}(t+\Delta t)\mathbf{U}(t)$ would trivially accomplish the goal of phase alignment, since then $\mathbf{U}(t+\Delta t)=\mathbf{U}(t)$. However, this $\bm{\Phi}^{\dagger}$ does not commute with $\mathbf{H}^{\text{diag}}(t+\Delta t)$. Hence, we propose the following ansatz for $\bm{\Phi}^{\dagger}$:

[TABLE]

Here, the operator $\hat{\mathcal{O}}$ ensures that $\bm{\Phi}^{\dagger}$ is orthogonal/unitary ($\bm{\Phi}^{\dagger}\bm{\Phi}=\mathbf{1}$) and $\hat{\mathcal{C}}$ ensures that $\bm{\Phi}^{\dagger}$ commutes with $\mathbf{H}^{\text{diag}}$ ($\mathbf{H}^{\text{diag}}\bm{\Phi}^{\dagger}=\bm{\Phi}^{\dagger}\mathbf{H}^{\text{diag}}$).

The commutation property can be achieved by defining

[TABLE]

where $\hat{\mathcal{P}}(\mathcal{S}_{i})$ is the projection operator on the $i$-th eigenspace $\mathcal{S}_{i}$ of $\mathbf{H}^{\text{diag}}(t+\Delta t)$:

[TABLE]

Here, the vectors $\mathbf{j}$ are the eigenvectors of $\mathbf{H}^{\text{diag}}(t+\Delta t)$. Because $\mathbf{H}^{\text{diag}}$ is by definition diagonal, the elements of $\hat{\mathcal{C}}\mathbf{V}$ can be written very compactly as

[TABLE]

It is easy to show that with this definition $\hat{\mathcal{C}}\mathbf{V}\mathbf{H}^{\text{diag}}=\mathbf{H}^{\text{diag}}\hat{\mathcal{C}}\mathbf{V}$.

The property $\bm{\Phi}^{\dagger}\bm{\Phi}=\mathbf{1}$ can be achieved if the operator $\hat{\mathcal{O}}$ orthonormalizes $\hat{\mathcal{C}}\mathbf{V}$. The optimal orthonormalization scheme here is Löwdin’s symmetric orthonormalization,203 since it changes $\hat{\mathcal{C}}\mathbf{V}$ as little as possible.204

The above-described algorithm leads for $\Delta t\rightarrow 0$ to $\mathbf{U}(t+\Delta t)\rightarrow\mathbf{U}(t)$. However, for finite $\Delta t$, one has to take some care. For the special case that the eigenvalues are non-degenerate ($H_{\alpha\alpha}^{\text{diag}}\neq H_{\beta\beta}^{\text{diag}}$ for all $\alpha,\beta$), one can show that

[TABLE]

Using $\tilde{U}_{\beta\alpha}(t+\Delta t)=U_{\beta\alpha}(t+\Delta t)\text{e}^{\text{i}\phi_{\alpha}}$ and reducing the $\delta_{i\alpha}$ gives:

[TABLE]

where $\mathbf{u}_{\alpha}$ is the $\alpha$-th column vector of $\mathbf{U}$. It is important to stress here that this only holds if $\mathbf{u}_{\alpha}(t+\Delta t)\cdot\mathbf{u}_{\alpha}(t)$ is real-valued, which will in general only be approximately true if no root flipping occurred between $t$ and $t+\Delta t$ (i.e., the $\alpha$-th eigenvector at time $t+\Delta t$ is the same as the $\alpha$-th eigenvector at time $t$). Thus, for finite timesteps $\Delta t$ we have to require that any root flipping is removed from $\mathbf{V}=\tilde{\mathbf{U}}^{\dagger}(t+\Delta t)\mathbf{U}(t)$. This can be accomplished by a simple algorithm that takes the matrix $\mathbf{V}$ and brings it into diagonally dominant form by reordering its columns.

Finally, it is crucial for the algorithm that $\Delta t$ is small. A criterion for smallness can be found from the antihermicity property of $\mathbf{U}^{\dagger}\frac{\mathrm{d}}{\mathrm{d}t}\mathbf{U}$. This term can be calculated by finite differences as

[TABLE]

The antihermicity of $\mathbf{U}^{\dagger}\frac{\mathrm{d}}{\mathrm{d}t}\mathbf{U}$ requires that its diagonal elements be zero, and hence it can be followed that the diagonal elements of $\mathbf{U}^{\dagger}(t+\Delta t)\mathbf{U}(t)$ need to be close to 1. Hence, we can choose a threshold $\varepsilon$ and require that

[TABLE]

If the condition is not fulfilled, a smaller value of $\Delta t$ is required.

The complete algorithm can thus be given by the following steps:

Numerically diagonalize: $\mathbf{H}^{\text{diag}}(t+\Delta t)=\tilde{\mathbf{U}}^{\dagger}(t+\Delta t)\mathbf{H}^{\text{MCH}}(t+\Delta t)\tilde{\mathbf{U}}(t+\Delta t)$. 2. 2.

Calculate $\mathbf{V}=\tilde{\mathbf{U}}^{\dagger}(t+\Delta t)\mathbf{U}(t)$. 3. 3.

Bring $\mathbf{V}$ into diagonally dominant form by exchanging columns. 4. 4.

Obtain the matrix $\hat{\mathcal{C}}\mathbf{V}$ from equation (49). 5. 5.

Orthonormalize $\hat{\mathcal{C}}\mathbf{V}$ by Löwdin’s symmetric orthonormalization.203 The result is $\bm{\Phi}^{\dagger}$. 6. 6.

Obtain the phase-controlled matrix $\mathbf{U}(t+\Delta t)=\tilde{\mathbf{U}}(t+\Delta t)\bm{\Phi}^{\dagger}$. 7. 7.

Check whether condition (53) holds. If not, divide the interval $[t,t+\Delta t]$ into subintervals, linearly interpolate $\mathbf{H}^{\text{MCH}}$ and carry out steps 1 to 6 for each subinterval. If the subintervals would be shorter than a given threshold, do not subdivide (termination of the recursive algorithm).

Finally, it should be stated that the algorithm outlined above also considers that the Hamiltonian can have degenerate eigenvalues. In this case, the mixing angles between the degenerate states are also not uniquely defined, in addition to the eigenvector phases. In step 4, the operator $\hat{\mathcal{C}}$ (see equation (49)) keeps off-diagonal elements of $\mathbf{V}$ corresponding to degenerate eigenfunctions, which leads to the removal of the arbitrariness in the mixing angles.