The change in electronic state from one spin multiplicity to another, known as intersystem crossing, occurs in molecules via the relativistic phenomenon of spin-orbit coupling. Current means of estimating intersystem crossing rates rely on the perturbative evaluation of spin-orbit coupling effects. This perturbative approach, valid in lighter atoms where spin-orbit coupling is weaker, is expected to break down for heavier elements where relativistic effects become dominant. Methods which incorporate spin-orbit effects variationally, such as the exact-two-component (X2C) method, will be necessary to treat this strong-coupling regime. We present a novel procedure which produces a diabatic basis of spin-pure electronic states coupled by spin-orbit terms, generated from fully variational relativistic calculations. This method is implemented within X2C using time-dependent density-functional theory and is compared to results from a perturbative relativistic study in the weak spin-orbit coupling regime. Additional calculations on a more strongly spin-orbit-coupled [UO2Cl4]2 complex further illustrate the strengths of this method. This procedure will be valuable in the estimation of intersystem crossing rates within strongly spin-coupled species.

Intersystem crossing (ISC) is the change in quantum state from one multiplicity to another. A formally spin-forbidden process, ISC nevertheless plays an important role in photochemistry.1 Due to exchange coupling, there are nearly always triplet excited states lower in energy than the lowest singlet excited state, and these states are by their spin state slow in decaying radiatively to the ground state. In systems where ISC is favorable, these low-lying and long-lasting triplet states can be efficiently populated, with direct consequences for the fluorescent and phosphorescent properties of a material. Accurate calculation of the ISC rates in a molecule, then, can be essential for understanding its photophysical properties.

With respect to spin eigenstates, intersystem crossing between different spin multiplicities occurs due to spin-orbit coupling (SOC), a fundamentally relativistic phenomenon.2 In addition to allowing for spin-forbidden transitions, relativistic effects include the contraction of core atomic orbitals in heavier elements, the splitting of degenerate multiplets, and many of the novel characteristics of transition metals, rare earth, and heavy elements.3,4 Modern electronic structure methods that include the important effects of spin-dependent quantum mechanics in chemical dynamics tend toward one of the two assumptions: they either treat spin coupling operators perturbatively via the spin-Hamiltonian or include the spin operators variationally via relativistic methods. While perturbative methods have their merits,5,6 particularly in the interpretation of ISC, they are limited in their description of the time-dependent nature of the spin-dependent many-electron wavefunction. This is especially true for phenomena driven by strong spin couplings, such as in heavy-metal complexes. A proper description of spin dynamics must arise from application of the time-dependent variational principle. However, within variational relativistic theory, interpretation and computation of ISC become challenging because spin is no longer a good quantum number. In this work, we present a novel procedure for calculating SOC matrix elements, which are needed for estimating ISC rates, within the framework of variational relativistic electronic structure theory.

To address intersystem crossing rates, it is necessary first to describe how systems may change quantum states. Within the Born-Oppenheimer approximation, the time-dependent electronic wavefunction Ψt,R(t) is a function of both time t and the time-dependent nuclear degrees of freedom R(t). This wavefunction may be expressed as

(1)

The set of electronic basis states ψI(R) will parametrically depend on the time-dependent R(t) but are otherwise assumed to be time-independent, and therefore, the time dependence of Ψ is located within the coefficients cI(t). The R(t) dependence of the basis states ψI(R) will be assumed for the remainder of the text. After inserting Eq. (1) into the time-dependent electronic Schrödinger equation, the equations of motion for the expansion coefficients are given by

(2)

Here, σLK(t)=ψL|t|ψK is the nonadiabatic coupling between states |ψL⟩ and |ψK⟩, and HLK(t)=ψL|ĤelR(t)|ψK is the electronic Hamiltonian expressed in the same basis. Population transfer between states is driven by the off-diagonal elements of H and σ in Eq. (2).

To simplify the dynamics, studies frequently perform simulations in one of the two simplifying bases. In the adiabatic basis, the basis states |ψI⟩ are taken from the eigenstates of Ĥel(R) at each point R(t) such that HLK(t) = ϵKδLK is diagonal, and only the nonadiabatic coupling matrix σLK retains off-diagonal elements. In the other simplifying case, the basis states are rotated to a diabatic basis |ψK′⟩ so as to minimize the nonadiabatic coupling term σL′K′ ∼ 0, leaving then only a nondiagonal HL′K′. As is difficult to diabatize states over all geometries for any but the smallest systems, an approximate diabatic basis is often chosen according to a variety of methods; for a review, see Ref. 7. Once a basis has been chosen, what remains is to sample the possible nuclear configurations R, for example by propagating them along classical trajectories.

When studying intersystem crossing, one is interested in the population transfer between electronic states of different spin multiplicities, most often between singlets and triplets.1,2 However, within nonrelativistic quantum mechanics, the electronic Hamiltonian does not contain spin-coupling terms capable of driving spin dynamics. Consequently, the integration of spin degrees of freedom over different spin states is zero, causing the nonadiabatic coupling σLK to vanish and preventing the population of one spin state from another.

When the spin-orbit coupling is relatively weak compared to the energy difference between spin states, perturbation theory can be used to introduce the ISC driving force into the nonrelativistic Hamiltonian. Among the most common is the Breit-Pauli spin-orbit operator,8–13 derived from the Dirac Hamiltonian, given as

(3)

where si is the spin vector for each electron. The first term contains single-electron SOC effects, while the second is a two-electron operator and is commonly neglected,12,13 parameterized into the one-electron term,14 or approximated using an atomic mean-field approach.15–17 After adding such an effective SO operator to the electronic Hamiltonian, states which previously had well-defined spins will couple together and no longer belong to a single multiplicity. In general, the Hamiltonian with the Breit-Pauli operator is unbounded from below, and naive eigenstate calculations will suffer from variational collapse.8,9 As a result, approximate eigenstates of the nonrelativistic Hamiltonian are typically calculated first, with the Breit-Pauli operator then applied perturbatively. If the Hamiltonian is again diagonalized in the resultant spin-perturbed basis, a set of adiabatic states naturally result.18–20 If the states instead remain spin-pure, they are often considered to be a nearly diabatic basis. This is because, as mentioned above, nonadiabatic coupling between spin-pure electronic states of different multiplicities is automatically zero. In contrast, nonadiabatic coupling between states of the same multiplicity, i.e. singlets, will not necessarily be zero. When singlet-triplet conversions are of particular interest, a first-order estimate of the ISC rate may be obtained by neglecting coupling between states of the same multiplicity.21–23 However, a full description of the intersystem crossing process will require the inclusion of all couplings: for example, efficient triplet-triplet internal conversion can increase the rate of ISC by preventing the reverse triplet-singlet conversion.

This approach, equivalent to a perturbation theory that is first order in spin-orbit coupling, is reasonable when treating lighter elements where SOC is weaker. Further down the periodic table, however, when SOC effects become increasingly strong, it is expected that approaches beyond first-order perturbation theory are necessary. In such cases, variational inclusion of relativistic effects is desirable and new diabatization techniques are needed for studies of ISC dynamics.

Within the kinetic-balance condition, we apply a spin separation technique so that the Dirac Hamiltonian can be written as8,9

(4)

where V is the scalar potential, c is the speed of light, and m is the electron mass. Spin and orbital angular momenta are coupled through the σ · pV × p term: the vector σ contains the Pauli spin matrices, and p is the linear momentum operator. Variational solutions to this Hamiltonian will no longer have a well-defined spin multiplicity, and the eigenstates will have four components. The first term in Eq. (4) is the spin-free portion of the Dirac Hamiltonian, which contains all scalar relativistic effects, while the second term gives rise to spin-couplings.

From the point of view of perturbation theory, the associated spin-diabatic states are those at zeroth order in terms of spin-coupling perturbation. That is, by using only the spin-free part of the Hamiltonian, variationally solving the Dirac equation will generate states that include scalar relativistic effects at the atomic orbital (AO) level but are still spin-pure, producing true singlets, triplets, etc. In the following discussion, {ψL} refer to variational solutions of the Dirac equation including all spin-couplings and {ψ̃K} are those obtained using only the spin-free portion of the Hamiltonian.

Unlike the perturbative approaches outlined above, the states {ψL} and {ψ̃K} are not simply related by the first-order spin-orbit coupling term, as such couplings are included in the infinite perturbation limit within the variational treatment. In order to calculate the couplings between states of different multiplicities, therefore, we search for a unitary transformation relating {ψL} and {ψ̃K}. The wavefunction overlap between the two sets of states is

(5)

which is non-Hermitian because {ψL} and {ψ̃K} are mutually nonorthogonal. To allow for a unitary transformation between the two sets, we perform a Löwdin orthogonalization of the spin-coupled states {ψL} within the basis of spin-free states {ψ̃K}. That is, we find the matrix U that diagonalizes the Hermitian matrix OO such that

(6)

where o is a diagonal matrix. The transformation T,24 

(7)

is subsequently a unitary transformation of the spin-coupled states, projecting them into the spin-free basis. Finally, taking the diagonal spin-coupled Hamiltonian H,

(8)

and rotating it into the spin-free basis

(9)

yields a diabatic Hamiltonian with off-diagonal couplings, whose eigenvalues are the same as those of spin-coupled H. In effect, this is a means of spin-purifying spin-orbit-coupled states.

Thus far, this procedure is entirely general for any method that can generate states from a variational relativistic electronic structure theory. However, the means of evaluating the wavefunction overlaps in Eq. (5) will be method-dependent. In this work, we will use the relativistic two-component time-dependent density-functional theory (TDDFT)25–27 method to illustrate the application of this approach to obtain the ISC coupling strength; for a recent review, see Ref. 28.

Relativistic electronic structure methods based on the Dirac equation employ a four-component wavefunction ansatz.8,9 For most chemically relevant studies, approximate two-component methods are usually of sufficient accuracy. The Douglas-Kroll-Hess (DKH)29–33 and exact-two-component (X2C)17,25,34–44 methods are popular approaches that decouple the four-component Hamiltonian and yield a reduced-dimension electronic two-component Hamiltonian, though they differ in the approximation used to decouple the large and small components of the relativistic wavefunction. Once an approximate decoupling is achieved, the two-component relativistic electronic Hamiltonian can be used with most electronic structure methods, provided a noncollinear spinor formalism is employed.45 It is of particular note that, unlike for effective SO operators like Breit-Pauli, X2C and DKH are bounded from below and do not suffer from variational collapse.8,9 Consequently, these two-component methods allow for the variational inclusion of relativistic effects at each stage of the calculation.

In this work, we use a direct atomic-orbital based X2C transformation with the torque-free spin-density approach.25,27,42 This transformation is only performed for the one-electron portion of the electronic Hamiltonian. To account partially for two-electron relativistic effects, the SO operator within X2C is scaled by a semiempirical fudge factor,14 reminiscent of the use of effective nuclear charges to account for nuclear screening of valence electrons by core electrons.

In the case of linear-response TDDFT, excited-state wavefunctions are not directly produced, but rather transition density matrices, complicating the evaluation of properties between excited states such as the wavefunction overlap in Eq. (5). Pseudowavefunctions, meant to reproduce approximately the linear-response eigenvectors, have been proposed,46,47 particularly in the context of evaluating nonadiabatic coupling vectors between excited states. For the latter method, it was shown that pseudowavefunctions produce coupling vectors in close agreement with those obtained from full quadratic response.48 If the Tamm-Dancoff approximation (TDA)49 is invoked [equivalent to configuration interaction singles (CIS) in the case of time-dependent Hartree-Fock], the auxiliary wavefunctions associated with excited states are expressed as linear combinations of Slater determinants singly excited from the DFT ground state,

(10)

where |Φ0⟩ is the ground-state determinant, and XiaK is the amplitude corresponding to the excitation from occupied orbital i to virtual orbital a for the state ψK. We invoke the TDA for two reasons. First, despite being slightly less accurate overall than full linear-response TDDFT, the TDA in nonrelativistic calculations has been shown to yield more accurate singlet-triplet energy gaps,50–52 which are vital for an accurate description of ISC. Second, in the absence of the so-called de-excitation amplitudes present in full linear-response TDDFT, the auxiliary wavefunctions are more intuitively constructed and are orthonormal, which greatly simplifies the procedure outlined here. In the case of TDA-TDDFT, the wavefunction overlap in Eq. (5) is therefore evaluated as

(11)

The expectation value Φ̃0|ãiãaâbâj|Φ0 is the overlap between two Slater determinants expressed in nonorthogonal sets of MOs. Such overlaps are commonly found both within nonorthogonal configuration interaction (NOCI)53 and whenever evaluating overlaps between wavefunctions at different time steps or nuclear geometries,24 and several methods for evaluating them efficiently have been discussed.54–56 The specific method employed in this work is outlined in the  Appendix.

It should be noted that a previous study57 evaluated the approximate overlap between excited states calculated using TDDFT at the scalar relativistic level and those with spin-orbit effects included at the zero-order regular approximation (ZORA)58 level of theory. This was done in order to estimate the composition of spin-orbit-coupled states in terms of spin-pure states, and in principle, that technique would also be capable of estimating SOC matrix elements in a similar fashion. However, as the previous work was principally interested in a qualitative comparison, it neglected rotations between occupied and virtual orbitals (in effect assuming that the two ground states were identical), an assumption justified only when SOC is very weak. In addition, X2C treats relativistic effects substantially more accurately than does ZORA, and as such, the method introduced here should yield more accurate couplings.

All calculations were performed using a locally modified copy of the developmental version of GAUSSIAN16.59 The geometry of the [Ru(bpy)3]2+ complex was optimized using the PBE functional,60 with the Ahlrichs double-ζ basis set with polarization on light atoms,61 and the LanL2DZ effective core potential and associated double-ζ basis set62 on Ru. X2C calculations used the same functional and basis set for the light atoms, with the relativistic all-electron Sapporo triple-ζ basis63 on Ru. The geometry of the [UO2Cl4]2 complex was optimized using the CAM-B3LYP functional,64 with the 6-31g(d) basis set on light atoms and the LanL2DZ effective core potential and associated double-ζ basis set on U. X2C calculations then used the cc-pVTZ65,66 basis set for light atoms and the relativistic SARC-DKH2 basis set67 for U. All nonrelativistic density functional theory calculations were performed on a pruned (75,302) grid (grid = UltraFine in GAUSSIAN16), while the convergence of relativistic calculations was improved by switching to a pruned (175,974) grid for light atoms and an unpruned (250,974) grid for heavy atoms (grid = SuperFine).

To test the validity of the approach outlined above for generating spin-orbit-coupling matrix elements, we compare our method against a recent DFT study68 of intersystem crossing in a ruthenium trisbipyridine ([Ru(bpy)3]2+) complex (Fig. 1), shown experimentally to have very rapid ISC.69,70 In the previous work, SOC matrix elements were evaluated by perturbatively applying the spin-orbit operator within the zero-order regular approximation (ZORA)58 to spin-pure states that included scalar relativistic effects at the ZORA level of theory.

FIG. 1.

Ruthenium trisbipyridine [Ru(bpy)3]2+ complex.

FIG. 1.

Ruthenium trisbipyridine [Ru(bpy)3]2+ complex.

Close modal

The energy levels of [Ru(bpy)3]2+, calculated using X2C, are plotted in Fig. 2. Due to the high symmetry of the molecule, there are several pairs of degenerate states in the spin-pure calculation. Additionally, within each set of triplets, the states are degenerate to within a millielectronvolt. In contrast, upon the inclusion of spin-orbit effects, there is a dramatic change in the energetics. Total energies shift by nearly a tenth of an electronvolt, and the degeneracy of some triplet states is broken by up to 30 meV. In this system, the relative energetic positions of the states make it easy to identify their predominant spin state; in general, however, once relativistic effects are included, it is no longer necessarily meaningful to speak of states as singlets or triplets, but rather they are linear superpositions of several multiplicities. Finally, Fig. 2 illustrates the results of using the diabatization scheme presented in this work. By rotating the spin-coupled states into the spin-pure basis, the degeneracy of triplet states is almost completely recovered: T4 has a splitting of 6 meV, while no other set of triplet states is split by more than 3 meV. Importantly, though, their relative energies are distinct from those of the initial spin-pure calculation; the states retain information about spin-orbit coupling effects, and as shown below, the resultant off-diagonal Hamiltonian elements contain the necessary information to evaluate intersystem crossing rates.

FIG. 2.

Energy level diagram of [Ru(bpy)3]2+, evaluated using X2C-TDDFT without spin-orbit coupling (left) and with spin-orbit coupling (middle). Spin-orbit effects shift states in energy and split the degeneracy of triplet states. Rotating spin-coupled states into the spin-pure basis (right) recovers much of the degeneracy of triplet states.

FIG. 2.

Energy level diagram of [Ru(bpy)3]2+, evaluated using X2C-TDDFT without spin-orbit coupling (left) and with spin-orbit coupling (middle). Spin-orbit effects shift states in energy and split the degeneracy of triplet states. Rotating spin-coupled states into the spin-pure basis (right) recovers much of the degeneracy of triplet states.

Close modal

The SOC matrix elements for the six lowest-energy singlets and nine lowest-energy triplets are compared in Table I: for degenerate pairs of states such as S2 and S3, the root-mean-square of both matrix elements is reported. Despite the SOC matrix elements spanning three orders of magnitude, we see close agreement between the two approaches. The mean absolute deviation for each singlet (or pair of singlets) is only a few tens of wavenumbers; over all 24 matrix elements, only one has a deviation greater than 80 cm−1. Given that ZORA uses different functional forms for both scalar and relativistic effects than X2C and also that SOC is evaluated perturbatively rather than variationally, it is not obvious that the two methods should yield similar values, yet the predictions for SOC in [Ru(bpy)3]2+ are in very close agreement. It appears likely, then, that this procedure is valid in the calculation of SOC matrix elements between singlets and triplets.

TABLE I.

Spin-orbit coupling matrix elements for [Ru(bpy)3]2+, evaluated in wavenumbers (cm−1) for the six lowest-energy singlets and nine lowest-energy triplets. For each singlet state, the triplets with the smallest and largest couplings are given, calculated using X2C (this work) and compared to those computed using first-order spin-orbit perturbed ZORA states.68 The mean absolute deviation (MAD) between the two methods is also given. For the complete data set, see the supplementary material.

Weakest couplingStrongest coupling
SingletTripletX2CZORATripletX2CZORAMAD
S0 T1 T2/T3 367 380 
S1 T1 T5/T6 626 626 15 
S2/S3 T1 22 29 T7/T8 632 622 23 
S4/S5 T4 17 23 T1 572 583 16 
Weakest couplingStrongest coupling
SingletTripletX2CZORATripletX2CZORAMAD
S0 T1 T2/T3 367 380 
S1 T1 T5/T6 626 626 15 
S2/S3 T1 22 29 T7/T8 632 622 23 
S4/S5 T4 17 23 T1 572 583 16 

To illustrate the effects of this method when applied to a more strongly spin-orbit-coupled system where perturbational methods might be expected to break down, we extend this analysis to the uranyl tetrachloride anion, [UO2Cl4]2. The actinide center leads to very strong SOC effects in this complex, and previous theoretical studies of this molecule exist for comparison.71,72 Excitation energies from spin-orbit-free calculations are given in Table II. When compared to results from a multireference CASPT2 calculation,71 very close agreement is seen, with most states agreeing to within a few tens of millielectronvolts. More striking are the results with spin-orbit coupling included, shown in Table III. Perturbative application of SOC in the CASPT2 calculation drops the energy of the lowest-lying excited states by less than 0.2 eV, while TDDFT with variational SOC whether two-component or four-component72 presents a massive change of nearly 0.5 eV. Clearly, variational treatment of SOC has a much larger effect in this complex than in the ruthenium complex above.

TABLE II.

Spin-orbit-free excitation energies of the six lowest-energy states in [UO2Cl4]2-, evaluated using X2C-TDDFT (this work) and CASPT2.71 Close agreement is seen between the two methods.

Excitation energy (eV)
StateX2C-TDDFTCASPT2
a3B1g 2.81 2.77 
a3B2g 2.90 2.87 
a3Eg 3.13 3.28 
a1Eg 3.58 3.71 
a1B1g 3.82 3.80 
a1B2g 3.89 3.90 
Excitation energy (eV)
StateX2C-TDDFTCASPT2
a3B1g 2.81 2.77 
a3B2g 2.90 2.87 
a3Eg 3.13 3.28 
a1Eg 3.58 3.71 
a1B1g 3.82 3.80 
a1B2g 3.89 3.90 
TABLE III.

Spin-orbit-coupled excitation energies of the 12 lowest-energy states in [UO2Cl4]2-, evaluated using X2C-TDDFT (this work), four-component TDDFT,72 and CASPT2.71 

Excitation energy (eV)
StateX2C-TDDFT4c-TDDFTCASPT2
a Eg 2.36 2.38 2.61 
a B2g 2.41 2.48 2.64 
a B1g 2.36 2.41 2.74 
b Eg 2.54 2.62 2.83 
b B1g 2.67 2.70 2.98 
b B2g 2.67 2.71 3.02 
c Eg 3.16 3.19 3.41 
a A2g 3.32 3.38 3.70 
a A1g 3.33 3.38 3.70 
d Eg 3.59 3.69 3.96 
c B1g 3.79 3.79 3.97 
c B2g 3.82 3.82 4.03 
Excitation energy (eV)
StateX2C-TDDFT4c-TDDFTCASPT2
a Eg 2.36 2.38 2.61 
a B2g 2.41 2.48 2.64 
a B1g 2.36 2.41 2.74 
b Eg 2.54 2.62 2.83 
b B1g 2.67 2.70 2.98 
b B2g 2.67 2.71 3.02 
c Eg 3.16 3.19 3.41 
a A2g 3.32 3.38 3.70 
a A1g 3.33 3.38 3.70 
d Eg 3.59 3.69 3.96 
c B1g 3.79 3.79 3.97 
c B2g 3.82 3.82 4.03 

The strengths of this method may be more clearly seen in the energy-level diagram of the [UO2Cl4]2 anion, shown in Fig. 3. Spin-orbit coupling leads to an enormous energy splitting of nearly half an eV, and the identification of states as singlets or triplets is clearly impossible. However, despite very powerful spin-orbit splitting, this method still recovers near-degeneracy of triplet states, with correspondingly large SOC matrix elements preserving the splitting of eigenvalues. The question remains of why exact degeneracy is not recovered. One possibility is that more excited states must be included in the sample space in order to adequately express spin-pure states in the spin-orbit-coupled basis, which is explored in Table IV. When only the first 12 states (primarily arising from the four lowest-energy triplets) are included, degeneracy is quite poor, with each given set of triplets possessing a standard deviation in energy of over 50 meV. However, including even the four next lowest-energy states (corresponding to the four lowest-energy singlets) dramatically improves the triplet degeneracies, and by including 35 states, degeneracy is recovered to the millielectronvolt level. Full degeneracy is still not completely recovered, demonstrating the incomplete nature of this transformation, but sufficient degeneracy is recovered that we believe this method will be applicable even to strongly coupled systems.

FIG. 3.

Energy level diagram of [UO2Cl4]2-, evaluated using X2C-TDDFT without spin-orbit coupling (left) and with spin-orbit coupling (middle). Despite energy splittings of nearly 0.5 eV upon the inclusion of SOC, near-degeneracy of triplet states is still recovered after rotating to a spin-pure basis (right). Dashed lines linking states are omitted for clarity due to very strong mixing between states.

FIG. 3.

Energy level diagram of [UO2Cl4]2-, evaluated using X2C-TDDFT without spin-orbit coupling (left) and with spin-orbit coupling (middle). Despite energy splittings of nearly 0.5 eV upon the inclusion of SOC, near-degeneracy of triplet states is still recovered after rotating to a spin-pure basis (right). Dashed lines linking states are omitted for clarity due to very strong mixing between states.

Close modal
TABLE IV.

Degeneracy of the three lowest spin-orbit-coupled triplets after being rotated into a spin-pure basis, as a function of the number of spin-pure states n included in the basis. The degeneracy is measured as the standard deviation σE of the energy of each set of triplets. More degeneracy is recovered as more states are added.

σE (eV)
Staten = 12n = 16n = 35
a3B1g 0.050 0.010 0.010 
a3B2g 0.056 0.007 0.007 
a3Eg 0.085 0.010 0.004 
σE (eV)
Staten = 12n = 16n = 35
a3B1g 0.050 0.010 0.010 
a3B2g 0.056 0.007 0.007 
a3Eg 0.085 0.010 0.004 

In this work, we have taken an important step toward the estimation of ISC rates using variational relativistic theory. While some perturbative treatments of spin-orbit coupling using scalar-relativistic methods have been proposed, for some of which18,73,74 analytical nuclear gradients have been implemented, the true strength of relativistic methods is their ability to include spin-orbit coupling effects fully variationally and at the molecular orbital level. However, variational treatment of ISC has not previously been implemented because (1) in the absence of a perturbative SO operator, it is nontrivial to evaluate couplings in a “spin-diabatic” basis and (2) the analytical excited-state nuclear gradients needed to evaluate nonadiabatic coupling75 are not yet available for two-component methods. The procedure outlined here, in contrast, takes states with spin-orbit effects included variationally and rotates them into a spin-pure basis. This effective spin-purification produces states of known multiplicities, making their identification easier, while also coupling them together in a manner that allows for intersystem crossing.

It should be noted that, in a finite atomic-centered basis and including only a subset of all possible excitations, the spin-pure states do not form a complete basis for the spin-coupled states, and care must be taken to ensure that not too much information is discarded in the change in the basis. In principle, the completeness of the two bases could be improved by increasing both the number of atomic basis functions and the number of excited states calculated, as shown in the uranyl tetrachloride example above. However, if the calculation is restricted to excited states with a finite number of excitations (e.g., singles, singles and doubles, etc.), the two sets of states may still not fully span one another, as they are comprised of different molecular orbitals. In practice, however, we have found that the ground state determinants overlap by at least 0.85 (0.99 in the case of the Ru complex above), allowing for the use of this method.

While SOC matrix elements are necessary in the treatment of intersystem crossing, more work is needed before variational relativistic theory can be used to estimate ISC rates. Needed next is the inclusion of nuclear dynamics effects. At present, in the absence of excited-state nuclear gradients, dynamics will need to be performed using a lower level of theory, a possibility explored in an upcoming publication. The development of such nuclear gradients would be invaluable to several extensions of relativistic theory: geometry optimizations, nuclear dynamics of spin-coupled electronic states, and of course the explicit evaluation of nonadiabatic coupling. Even without such gradients, however, the machinery for calculating wavefunction overlaps in the manner introduced here could be used to probe coupling in the adiabatic basis by the numerical evaluation of ψK|t|ψL,76,77 or by implementing the local diabatization method,24,78 where states at one time step are rotated into the basis of states at the preceding time step. Also to be determined is whether spin-pure states must be evaluated at each time step of a trajectory or whether a single representative set (or at least a smaller number of sets) of spin-pure states could provide a sufficiently accurate diabatic basis. By combining the evaluation of spin-orbit coupling using variational relativistic theory with a method for incorporating nuclear dynamics, we will be able to directly model population transfer between states of different multiplicities after photoexcitation, allowing us to accurately and efficiently estimate intersystem crossing rates.

The full set of spin-orbit coupling matrix elements for the ruthenium complex calculated using X2C-TDDFT and their comparison to those from ZORA68 are presented in the supplementary material.

The development of the two-component electronic structure method is funded by the U.S. Department of Energy (Grant No. DE-SC0006863). The development of nonperturbative spectroscopic methods is supported by the National Science Foundation (Grant No. CHE-1856210). Computational intersystem crossing is partially supported by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences, Division of Chemical Sciences, Geosciences and Biosciences, through Argonne National Laboratory under Contract No. DE-AC02-06CH11357. Computations were facilitated through the use of advanced computational, storage, and networking infrastructure provided by the Hyak supercomputer system at the University of Washington, funded by the Student Technology Fee and the National Science Foundation (Grant No. MRI-1624430).

For two sets of nonorthogonal molecular orbitals (MOs) defined by the matrices C̃ and C, constructed from the same atomic orbital (AO) basis, the overlap between the MOs |ϕ̃p and |ϕq⟩ is given by

(A1)

where S is the AO overlap matrix. The wavefunction overlap between two n-electron Slater determinants |Φ̃ and |Φ⟩ constructed from these MOs is given by the determinant

(A2)

Naively implemented, the overlaps in Eq. (11) would require calculating the overlap Φ̃ia|Φjb for all n2v2 possibilities i,j,a,b, where |Φia is the Slater determinant formed by promoting an electron from occupied orbital i to virtual orbital a. However, a substantial reduction in complexity can be realized by noting that rotations among occupied orbitals do not change the energy of the DFT ground state. One is free, then, to rotate the occupied orbitals C̃O and CO, respectively, by the unitary transformations U and V, found by performing a singular value decomposition (SVD) of the ground-state overlap

(A3)

where s is a diagonal matrix. The overlap between the respective ground state determinants Φ̃0|Φ0 is now (to within a phase) equal to the product of the singular values ∏si. The new occupied orbitals, given as C̃O=C̃OU and CO = COV, in the two sets are then biorthonormal: that is, ϕĩ|ϕj=δij. By rotating M into the new basis,

(A4)

as well as the TDA-TDDFT excitation coefficients,

(A5)

the overlap Φ̃ia|Φjb can now much more simply be evaluated as

(A6)

reducing the complexity of the overall overlap calculation from O(n5v2) to O(n2v2).

1.
C. M.
Marian
, “
Spin-orbit coupling and intersystem crossing in molecules
,”
Wiley Interdiscip. Rev.: Comput. Mol. Sci.
2
,
187
203
(
2012
).
2.
T. J.
Penfold
,
E.
Gindensperger
,
C.
Daniel
, and
C. M.
Marian
, “
Spin-vibronic mechanism for intersystem crossing
,”
Chem. Rev.
118
,
6975
7025
(
2018
).
3.
P.
Pyykkö
, “
Relativistic effects in structural chemistry
,”
Chem. Rev.
88
,
563
294
(
1988
).
4.
P.
Pyykkö
, “
Relativistic effects in chemistry: More common than you thought
,”
Annu. Rev. Phys. Chem.
63
,
45
64
(
2012
).
5.
F.
Neese
, “
Importance of direct spin-spin coupling and spin-flip excitations for the zero-field splittings of transition metal complexes: A case study
,”
J. Am. Chem. Soc.
128
,
10213
10222
(
2006
).
6.
M.
Seth
and
T.
Ziegler
, “
Formulation of magnetically perturbed time-dependent density functional theory
,”
J. Chem. Phys.
127
,
134108
(
2007
).
7.
T.
Van Voorhis
,
T.
Kowalczyk
,
B.
Kaduk
,
L.-P.
Wang
,
C.-L.
Cheng
, and
Q.
Wu
, “
The diabatic picture of electron transfer, reaction barriers, and molecular dynamics
,”
Annu. Rev. Phys. Chem.
61
,
149
170
(
2010
).
8.
K. G.
Dyall
and
K.
Fægri
, Jr.
,
Introduction to Relativistic Quantum Chemistry
(
Oxford University Press
,
2007
).
9.
M.
Reiher
and
A.
Wolf
,
Relativistic Quantum Chemistry
, 2nd ed. (
Wiley-VCH
,
2015
).
10.
S.
Chiodo
and
N.
Russo
, “
Determination of spin-orbit coupling contributions in the framework of density functional theory
,”
J. Comput. Chem.
29
,
912
920
(
2008
).
11.
S. G.
Chiodo
and
N.
Russo
, “
DFT spin-orbit coupling between singlet and triplet excited states: A case of psoralen compounds
,”
Chem. Phys. Lett.
490
,
90
96
(
2010
).
12.
D. G.
Fedorov
and
M. S.
Gordon
, “
A study of the relative importance of one and two-electron contributions to spin-orbit coupling
,”
J. Chem. Phys.
112
,
5611
5623
(
2000
).
13.
X.
Gao
,
S.
Bai
,
D.
Fazzi
,
T.
Niehaus
,
M.
Barbatti
, and
W.
Thiel
, “
Evaluation of spin-orbit couplings with linear-response time-dependent density functional methods
,”
J. Chem. Theory Comput.
13
,
515
524
(
2017
).
14.
J. C.
Boettger
, “
Approximate two-electron spin-orbit coupling term for density-functional-theory DFT calculations using the Douglas-Kroll-Hess transformation
,”
Phys. Rev. B
62
,
7809
7815
(
2000
).
15.
B. A.
Heß
,
C. M.
Marian
,
U.
Wahlgren
, and
O.
Gropen
, “
A mean-field spin-orbit method applicable to correlated wavefunctions
,”
J. Chem. Phys.
251
,
365
371
(
1996
).
16.
L.
Cheng
and
J.
Gauss
, “
Perturbative treatment of spin-orbit coupling within spin-free exact two-component theory
,”
J. Chem. Phys.
141
,
164107
(
2014
).
17.
Z.
Li
,
Y.
Xiao
, and
W.
Liu
, “
On the spin separation of algebraic two-component relativistic Hamiltonians
,”
J. Chem. Phys.
137
,
154114
(
2012
).
18.
N.
Bellonzi
,
G. R.
Medders
,
E.
Epifanovsky
, and
J. E.
Subotnik
, “
Configuration interaction singles with spin-orbit coupling: Constructing spin-adiabatic states and their analytical nuclear gradients
,”
J. Chem. Phys.
150
,
014106
(
2019
).
19.
S.
Mai
,
P.
Marquetand
, and
L.
González
, “
Nonadiabatic dynamics: The SHARC approach
,”
Wiley Interdiscip. Rev.: Comput. Mol. Sci.
8
,
e1370
(
2018
).
20.
S.
Mai
,
F.
Plasser
,
M.
Pabst
,
F.
Neese
,
A.
Köhn
, and
L.
González
, “
Surface hopping dynamics including intersystem crossing using the algebraic diagrammatic construction method
,”
J. Chem. Phys.
147
,
184109
(
2017
).
21.
G.
Granucci
,
M.
Persico
, and
G.
Spighi
, “
Surface hopping trajectory simulations with spin-orbit and dynamical couplings
,”
J. Chem. Phys.
137
,
22A501
(
2012
).
22.
F.
de Carvalho
and
I.
Tavernelli
, “
Nonadiabatic dynamics with intersystem crossings: A time-dependent density functional theory implementation
,”
J. Chem. Phys.
143
,
224105
(
2015
).
23.
G.
Cui
and
W.
Thiel
, “
Generalized trajectory surface-hopping method for internal conversion and intersystem crossing
,”
J. Chem. Phys.
141
,
124101
(
2014
).
24.
G.
Granucci
,
M.
Persico
, and
A.
Toniolo
, “
Direct semiclassical simulation of photochemical processes with semiempirical wave functions
,”
J. Chem. Phys.
114
,
10608
10615
(
2001
).
25.
F.
Egidi
,
S.
Sun
,
J. J.
Goings
,
G.
Scalmani
,
M. J.
Frisch
, and
X.
Li
, “
Two-component non-collinear time-dependent spin density functional theory for excited state calculations
,”
J. Chem. Theory Comput.
13
,
2591
2603
(
2017
).
26.
J. M.
Kasper
,
P. J.
Lestrange
,
T. F.
Stetina
, and
X.
Li
, “
Modeling l2,3-edge x-ray absorption spectroscopy with real-time exact two-component relativistic time-dependent density functional theory
,”
J. Chem. Theory Comput.
14
,
1998
2006
(
2018
).
27.
A.
Petrone
,
D. B.
Williams-Young
,
S.
Sun
,
T. F.
Stetina
, and
X.
Li
, “
An efficient implementation of two-component relativistic density functional theory with torque-free auxiliary variables
,”
Eur. Phys. J. B
91
,
169
(
2018
).
28.
W.
Liu
and
Y.
Xiao
, “
Relativistic time-dependent density functional theories
,”
Chem. Soc. Rev.
47
,
4481
4509
(
2018
).
29.
M.
Douglas
and
N. M.
Kroll
, “
Quantum electrodynamical corrections to the fine structure of helium
,”
Ann. Phys.
82
,
89
155
(
1974
).
30.
B. A.
Hess
, “
Relativistic electronic-structure calculations employing a two-component no-pair formalism with external-field projection operators
,”
Phys. Rev. A
33
,
3742
3748
(
1986
).
31.
M.
Reiher
and
A.
Wolf
, “
Exact decoupling of the Dirac Hamiltonian. I. General theory
,”
J. Chem. Phys.
121
,
2037
(
2004
).
32.
M.
Reiher
and
A.
Wolf
, “
Exact decoupling of the Dirac Hamiltonian. II. The generalized Douglas-Kroll-Hess transformation up to arbitrary order
,”
J. Chem. Phys.
121
,
10945
(
2004
).
33.
A.
Wolf
and
M.
Reiher
, “
Exact decoupling of the Dirac Hamiltonian. III. Molecular properties
,”
J. Chem. Phys.
121
,
064102
(
2006
).
34.
W.
Kutzelnigg
and
W.
Liu
, “
Quasirelativistic theory equivalent to fully relativistic theory
,”
J. Chem. Phys.
123
,
241102
(
2005
).
35.
W.
Liu
and
D.
Peng
, “
Infinite-order quasirelativistic density functional method based on the exact matrix quasirelativistic theory
,”
J. Chem. Phys.
125
,
044102
(
2006
).
36.
D.
Peng
,
W.
Liu
,
Y.
Xiao
, and
L.
Cheng
, “
Making four- and two-component relativistic density functional methods fully equivalent based on the idea of ‘from atoms to molecule
,’”
J. Chem. Phys.
127
,
104106
(
2007
).
37.
M.
Ilias
and
T.
Saue
, “
An infinite-order relativistic Hamiltonian by a simple one-step transformation
,”
J. Chem. Phys.
126
,
064102
(
2007
).
38.
W.
Liu
and
D.
Peng
, “
Exact two-component Hamiltonians revisited
,”
J. Chem. Phys.
131
,
031104
(
2009
).
39.
W.
Liu
, “
Ideas of relativistic quantum chemistry
,”
Mol. Phys.
108
,
1679
1706
(
2010
).
40.
T.
Saue
, “
Relativistic Hamiltonians for chemistry: A primer
,”
ChemPhysChem
12
,
3077
3094
(
2011
).
41.
D.
Peng
,
N.
Middendorf
,
F.
Weigend
, and
M.
Reiher
, “
An efficient implementation of two-component relativistic exact-decoupling methods for large molecules
,”
J. Chem. Phys.
138
,
184105
(
2013
).
42.
F.
Egidi
,
J. J.
Goings
,
M. J.
Frisch
, and
X.
Li
, “
Direct atomic-orbital-based relativistic two-component linear response method for calculating excited-state fine structures
,”
J. Chem. Theory Comput.
12
,
3711
3718
(
2016
).
43.
J. J.
Goings
,
J. M.
Kasper
,
F.
Egidi
,
S.
Sun
, and
X.
Li
, “
Real time propagation of the exact two component time-dependent density functional theory
,”
J. Chem. Phys.
145
,
104107
(
2016
).
44.
L.
Konecny
,
M.
Kadek
,
S.
Komorovsky
,
O. L.
Malkina
,
K.
Ruud
, and
M.
Repisky
, “
Acceleration of relativistic electron dynamics by means of X2C transformation: Application to the calculation of nonlinear optical properties
,”
J. Chem. Theory Comput.
12
,
5823
5833
(
2016
).
45.
J. J.
Goings
,
F.
Egidi
, and
X.
Li
, “
Current development of non-collinear electronic structure theory
,”
Int. J. Quantum Chem.
118
,
e25398
(
2018
).
46.
I.
Tavernelli
,
B. F. E.
Curchod
,
A.
Laktionov
, and
U.
Rothlisberger
, “
Nonadiabatic coupling vectors for excited states within time-dependent density functional theory in the Tamm-Dancoff approximation and beyond
,”
J. Chem. Phys.
133
,
194104
(
2010
).
47.
E. C.
Alguire
,
Q.
Ou
, and
J. E.
Subotnik
, “
Calculating derivative couplings between time-dependent Hartree-Fock excited states with pseudo-wavefunctions
,”
J. Phys. Chem. B
119
,
7140
7149
(
2015
).
48.
X.
Zhang
and
J. M.
Herbert
, “
Analytic derivative couplings in time-dependent density functional theory: Quadratic response theory versus pseudo-wavefunction approach
,”
J. Chem. Phys.
142
,
064109
(
2015
).
49.
S.
Hirata
and
M.
Head-Gordon
, “
Time-dependent density functional theory within the Tamm-Dancoff approximation
,”
Chem. Phys. Lett.
314
,
291
299
(
1999
).
50.
F.
Dinkelbach
,
M.
Kleinschmidt
, and
C. M.
Marian
, “
Assessment of interstate spin-orbit couplings from linear response amplitudes
,”
J. Chem. Theory Comput.
13
,
749
766
(
2017
).
51.
M. J. G.
Peach
,
M. J.
Williamson
, and
D. J.
Tozer
, “
Influence of triplet instabilities in TDDFT
,”
J. Chem. Theory Comput.
7
,
3578
3585
(
2011
).
52.
A. J.
Atkins
,
F.
Talotta
,
L.
Freitag
,
M.
Boggio-Pasqua
, and
L.
González
, “
Assessing excited state energy gaps with time-dependent density functional theory on Ru(II) complexes
,”
J. Chem. Theory Comput.
13
,
4123
4145
(
2017
).
53.
H.
Goto
,
M.
Kojo
,
A.
Sasaki
, and
K.
Hirose
, “
Essentially exact ground-state calculations by superpositions of nonorthogonal Slater determinants
,”
Nanoscale Res. Lett.
8
,
200
(
2013
).
54.
F.
Plasser
,
M.
Ruckenbauer
,
S.
Mai
,
M.
Oppel
,
P.
Marquetand
, and
L.
González
, “
Efficient and flexible computation of many-electron wave function overlaps
,”
J. Chem. Theory Comput.
12
,
1207
1219
(
2016
).
55.
M.
Sapunar
,
T.
Piteša
,
D.
Davidović
, and
N.
Došlić
, “
Highly efficient algorithms for CIS type excited state wave function overlaps
,”
J. Chem. Theory Comput.
15
,
3461
3469
(
2019
).
56.
G. A.
Meek
and
B. G.
Levine
, “
Evaluation of the time-derivative coupling for accurate electronic state transition probabilities from numerical simulations
,”
J. Phys. Chem. Lett.
5
,
2351
2356
(
2014
).
57.
F.
Wang
and
T.
Ziegler
, “
Theoretical study of the electronic spectra of square-planar platinum (II) complexes based on the two-component relativistic time-dependent density-functional theory
,”
J. Chem. Phys.
123
,
194102
(
2005
).
58.
E.
van Lenthe
,
E. J.
Baerends
, and
J. G.
Snijders
, “
Relativistic regular two-component Hamiltonians
,”
J. Chem. Phys.
99
,
4597
4610
(
1993
).
59.
M. J.
Frisch
,
G. W.
Trucks
,
H. B.
Schlegel
,
G. E.
Scuseria
,
M. A.
Robb
,
J. R.
Cheeseman
,
G.
Scalmani
,
V.
Barone
,
G. A.
Petersson
,
H.
Nakatsuji
,
X.
Li
,
M.
Caricato
,
A. V.
Marenich
,
J.
Bloino
,
B. G.
Janesko
,
R.
Gomperts
,
B.
Mennucci
,
H. P.
Hratchian
,
J. V.
Ortiz
,
A. F.
Izmaylov
,
J. L.
Sonnenberg
,
D.
Williams-Young
,
F.
Ding
,
F.
Lipparini
,
F.
Egidi
,
J.
Goings
,
B.
Peng
,
A.
Petrone
,
T.
Henderson
,
D.
Ranasinghe
,
V. G.
Zakrzewski
,
J.
Gao
,
N.
Rega
,
G.
Zheng
,
W.
Liang
,
M.
Hada
,
M.
Ehara
,
K.
Toyota
,
R.
Fukuda
,
J.
Hasegawa
,
M.
Ishida
,
T.
Nakajima
,
Y.
Honda
,
O.
Kitao
,
H.
Nakai
,
T.
Vreven
,
K.
Throssell
,
J. A.
Montgomery
, Jr.
,
J. E.
Peralta
,
F.
Ogliaro
,
M. J.
Bearpark
,
J. J.
Heyd
,
E. N.
Brothers
,
K. N.
Kudin
,
V. N.
Staroverov
,
T. A.
Keith
,
R.
Kobayashi
,
J.
Normand
,
K.
Raghavachari
,
A. P.
Rendell
,
J. C.
Burant
,
S. S.
Iyengar
,
J.
Tomasi
,
M.
Cossi
,
J. M.
Millam
,
M.
Klene
,
C.
Adamo
,
R.
Cammi
,
J. W.
Ochterski
,
R. L.
Martin
,
K.
Morokuma
,
O.
Farkas
,
J. B.
Foresman
, and
D. J.
Fox
, gaussian 16, Revision B.01,
Gaussian, Inc.
,
Wallingford, CT
,
2016
.
60.
J. P.
Perdew
,
K.
Burke
, and
M.
Ernzerhof
, “
Generalized gradient approximation made simple
,”
Phys. Rev. Lett.
77
,
3865
3868
(
1996
).
61.
F.
Weigend
and
R.
Ahlrichs
, “
Balanced basis sets of split valence, triple zeta valence and quadruple zeta valence quality for H to Rn: Design and assessment of accuracy
,”
Phys. Chem. Chem. Phys.
7
,
3297
(
2005
).
62.
P. J.
Hay
and
W. R.
Wadt
, “
Ab initio effective core potentials for molecular calculations. Potentials for K to Au including the outermost core orbitals
,”
J. Chem. Phys.
82
,
299
310
(
1985
).
63.
T.
Noro
,
M.
Sekiya
, and
T.
Koga
, “
Segmented contracted basis sets for atoms H through Xe: Sapporo-(DK)-nZP sets (n = D, T, Q)
,”
Theor. Chem. Acc.
131
,
1124
(
2012
).
64.
T.
Yanai
,
D. P.
Tew
, and
N. C.
Handy
, “
A new hybrid exchange-correlation functional using the coulomb-attenuating method (CAM-B3LYP)
,”
Chem. Phys. Lett.
393
,
51
57
(
2004
).
65.
T. H.
Dunning
, “
Gaussian basis sets for use in correlated molecular calculations. I. The atoms boron through neon and hydrogen
,”
J. Chem. Phys.
90
,
1007
1023
(
1989
).
66.
D. E.
Woon
and
T. H.
Dunning
, “
Gaussian basis sets for use in correlated molecular calculations. III. The atoms aluminum through argon
,”
J. Chem. Phys.
98
,
1358
1371
(
1993
).
67.
D. A.
Pantazis
and
F.
Neese
, “
All-electron scalar relativistic basis sets for the actinides
,”
J. Chem. Theory Comput.
7
,
677
684
(
2011
).
68.
A. J.
Atkins
and
L.
González
, “
Trajectory surface-hopping dynamics including intersystem crossing in [Ru(bpy)3]2+
,”
J. Phys. Chem. Lett.
8
,
3840
3845
(
2017
).
69.
A.
Cannizzo
,
F.
van Mourik
,
W.
Gawelda
,
G.
Zgrablic
,
C.
Bressler
, and
M.
Chergui
, “
Broadband femtosecond fluorescence spectroscopy of [Ru(bpy)3]2+.
,”
Angew. Chem., Int. Ed.
45
,
3174
3176
(
2006
).
70.
N. H.
Damrauer
,
G.
Cerullo
,
A.
Yeh
,
T. R.
Boussie
,
C. V.
Shank
, and
J. K.
McCusker
, “
Femtosecond dynamics of excited-state evolution in [Ru(bpy)3]2+.
,”
Science
275
,
54
57
(
1997
).
71.
K.
Pierloot
and
E.
van Besien
, “
Electronic spectrum and spectrum of UO  22+ and [UO2Cl4]2−
,”
J. Chem. Phys.
123
,
204309
(
2005
).
72.
P.
Tecmer
,
R.
Bast
,
K.
Ruud
, and
L.
Visscher
, “
Charge-transfer excitations in uranyl tetrachloride ([UO2Cl4]2): How reliable are electronic spectra from relativistic time-dependent density functional theory?
,”
J. Phys. Chem. A
116
,
7397
7404
(
2012
).
73.
X.
Zhou
,
Z.
Cao
, and
F.
Wang
, “
Analytical energy gradients for ionized states using equation-of-motion coupled-cluster theory with spin-orbit coupling
,”
J. Chem. Phys.
150
,
154114
(
2019
).
74.
G.
Granucci
and
M.
Persico
, “
Gradients for configuration interaction energies with spin-orbit coupling in a semiempirical framework
,”
J. Comput. Chem.
32
,
2690
2696
(
2011
).
75.
J. C.
Tully
, “
Perspective: Nonadiabatic dynamics theory
,”
J. Chem. Phys.
137
,
22A301
(
2012
).
76.
R.
Mitrić
,
U.
Werner
, and
V.
Bonačić-Koutecký
, “
Nonadiabatic dynamics and simulation of time resolved photoelectron spectra within time-dependent density functional theory: Ultrafast photoswitching in benzylideneaniline
,”
J. Chem. Phys.
129
,
164118
(
2008
).
77.
S.
Hammes-Schiffer
and
J. C.
Tully
, “
Proton transfer in solution: Molecular dynamics with quantum transitions
,”
J. Chem. Phys.
101
,
4657
4667
(
1994
).
78.
F.
Plasser
,
G.
Granucci
,
J.
Pittner
,
M.
Barbatti
,
M.
Persico
, and
H.
Lischka
, “
Surface hopping dynamics using a locally diabatic formalism: Charge transfer in the ethylene dimer cation and excited state dynamics in the 2-pyridone dimer
,”
J. Chem. Phys.
137
,
22A514
(
2012
).

Supplementary Material