We evaluate the effect of electronic decoherence on intersystem crossing in the photodynamics of thioformaldehyde. First, we show that the state-averaged complete-active-space self-consistent field electronic structure calculations with a properly chosen active space of 12 active electrons in 10 active orbitals can predict the potential energy surfaces and the singlet–triplet spin–orbit couplings quite well for CH2S, and we use this method for direct dynamics by coherent switching with decay of mixing (CSDM). We obtain similar dynamical results with CSDM or by adding energy-based decoherence to trajectory surface hopping, with the population of triplet states tending to a small steady-state value over 500 fs. Without decoherence, the state populations calculated by the conventional trajectory surface hopping method or the semiclassical Ehrenfest method gradually increase. This difference shows that decoherence changes the nature of the results not just quantitatively but qualitatively.

Intersystem crossing (ISC) is a nonradiative transition between different electronic spin states, and it is crucial for many applications, e.g., gas-phase reactions of unsaturated hydrocarbons with an oxygen atom in its ground 3P state1–4 and processes in energy-related materials.5–13 However, it is difficult to simulate. Population transfer between different spin states is due to spin–orbit coupling (SOC); it becomes very important when the mass of one or more of the involved atoms is increased, but it can also be important for atoms as light as sulfur, and for some small organic molecules, ISC happens rapidly—on a subpicosecond time scale.14–19 

Simulation of nonadiabatic processes involving ISC requires appropriate treatment of nonadiabatic transitions induced not only by SOC, which causes ISC, but also by spin-conserving nuclear-momentum coupling, which is responsible for non-Born–Oppenheimer internal conversion (IC). The Sharc program20–23 was designed to handle both ISC and IC in a consistent way by propagating semiclassical trajectories on a fully diagonalized electronic Hamiltonian, which is called the diagonal representation. The diagonal representation is constructed by diagonalization of the blocked molecular-Coulomb–Hamiltonian (MCH), for which the diagonal blocks are for the states of a given spin, and the off-diagonal blocks connect states with different spins. (Alternative labels24 for these representations are fully adiabatic for diagonal and valence adiabatic for MCH.)

The Sharc program was first developed for trajectory surface hopping (TSH),25–30 which replaces the nuclear wave packet by a swarm of independent trajectories that propagate between hops on single fully adiabatic potential energy surfaces. Nonradiative transitions are achieved in TSH by stochastic hopping of the trajectories from one surface to another; decoherence is neglected. Recently, we have extended the capability of Sharc to treat mixed quantum-classical dynamics with a self-consistent potential (SCP) and decoherence by using the coherent switching with decay of mixing (CSDM)31 method. We also implemented the semiclassical Ehrenfest (SE)32–35 method that employs the SCP without decoherence. The SCP is a mean field corresponding to a coherent superposition of several adiabatic electronic states, and such a method is sometimes called the mean-field method.36,37 Both the CSDM method and the SE method, like TSH, involve batches of independent semiclassical trajectories. The semiclassical trajectories involve classical propagation of nuclear motion, but they are not purely classical because they either involve hopping between two or more surfaces (in TSH) or propagation on a SCP that is derived from two or more quantal electronic states (in CSDM or SE).

The mean-field approximation underlying the SCP is appealing because it accounts for the coherence in strong-interaction regions, i.e., regions in the vicinity of locally avoided crossings, conical intersections, and crossings of different spin states. As a consequence, SCP methods are robust on the choice of representation (i.e., the statistical results from an ensemble of trajectories are almost the same for trajectories propagating on adiabatic or diabatic potential surfaces if one simulates IC processes, or on diagonal or MCH surfaces when one simulates ISC processes); this is consistent with fully quantal dynamics. The SE approximation suffers from overcoherence,41–42 which leads to unphysical trajectories for long-time dynamics. CSDM introduces decoherence into SE by adding a non-Markovian decay of the off-diagonal elements of the electronic density matrix. Between strong-interaction regions, the electronic density matrix of CSDM tends to become diagonal in the pointer basis on a time scale called the decoherence time. Algorithmically, when treating ISC, one may choose either the fully adiabatic basis or the MCH basis as the pointer basis. The CSDM is consistent with the Liouville–von Neumann equation of motion.43 

Standard TSH, like SE, suffers from overcoherence. One can introduce decoherence into TSH in various ways. For example, within the independent-trajectory scheme, one can employ an overlap-based decoherence correction (TSH-ODC),44 an energy-based decoherence correction (TSH-EDC)45,46 (based on an earlier approximation31 to the decoherence time in terms of energy gaps and nuclear kinetic energies), or the fewest-switches time-uncertainty stochastic-decoherence (FSTU/SD)47 extension of TSH.

In addition to requiring a choice of dynamical method, simulations of photophysical dynamics also require a choice of electronic structure method to compute potential energy surfaces and state couplings. Mai et al.23 recently examined the effect of choosing different electronic structure methods on the photodynamics of thioformaldehyde. They found quite different results with state-averaged complete-active-space self-consistent field48 (SA-CASSCF) calculations and multistate complete active space second-order perturbation theory49,50 (MS-CASPT2), and this motivates further examination of the adequacy of SA-CASSCF.

With the aim of investigating the important role of the electronic structure method and the inclusion of decoherence in intersystem crossing dynamics, we here report a theoretical and computational study of the photodynamics of thioformaldehyde (CH2S), which has been the subject of several spectroscopic investigations.51–53 We report simulations by CSDM, SE, TSH, and TSH-EDC. This is the first example of applying the CSDM and SE capability of Sharc to the study of an ISC process.

Second-order Møller–Plesset perturbation theory54 (MP2) and SA-CASSCF calculations were performed with the Molpro55,56 software package; MS-CASPT2 calculations carried out by Mai et al.23 were done using the Molcas57 software package. Orbital plots were made with the help of Multiwfn58 and Vmd59 software. Multireference configuration interactions including single and double excitations with the Pople size-extensivity correction60,61 (MRCISD+P) are taken from Ref. 23.

Semiclassical trajectory calculations were performed with a locally modified version36,37,62 of the Sharc software package,22,63 which employs direct dynamics, which in the present context denotes that electronic structure calculations are calculated presently as energies, gradients, and couplings are needed during the propagation of the trajectories.

The ground-state equilibrium structure of thioformaldehyde was optimized, and its frequencies were calculated by MP2 with the 6-31G*64,65 basis set. We will also present some results for geometries in which the C–S bond is stretched from its equilibrium length (RC–S = 1.617 Å) with the other internal coordinates fixed at their values in the equilibrium structure.

The ground and excited singlet and triplet states of the CH2S were first investigated with the SA-CASSCF with the 6-31G* basis set and active space of 12 electrons in 10 orbitals, denoted as SA-CASSCF(12,10), as shown in Fig. 1 and Fig. S1 of the supplementary material for the geometry with the C–S distance increased to RC–S = 2.10 Å. A small active space of 10 electrons in 6 orbitals (the encased orbitals in Fig. 1) was also tested, denoted as SA-CASSCF(10,6). The state-averaging in the SA-CASSCF calculations includes two singlet (S0, S1) and two triplet (T1, T2) states, denoted as SA(2S+2T)-CASSCF(12,10), and was performed with six frozen core orbitals and with the second-order Douglass–Kroll–Hess scalar relativistic Hamiltonian66,67 (by setting DHKO = 2 in Molpro).

FIG. 1.

The ten state-averaged natural active orbitals used in the CASSCF(12,10) calculations at the S0 equilibrium geometry. The encased orbitals (1 and 3–7) are the ones involved in the small active space CASSCF(10,6) calculations.

FIG. 1.

The ten state-averaged natural active orbitals used in the CASSCF(12,10) calculations at the S0 equilibrium geometry. The encased orbitals (1 and 3–7) are the ones involved in the small active space CASSCF(10,6) calculations.

Close modal

In Secs. III and IV, we abbreviate MP2/6-31G* as MP2, SA(2S+2T)-CASSCF(12,10)/6-31G* as CASSCF(12,10), and SA(2S+2T)-CASSCF(10,6)/6-31G* as CASSCF(10,6).

The magnitude of the S1–T2 spin–orbit coupling will be illustrated by a strength parameter defined by

SOC=M=1,0,+1S1|HSOC|T2M2,
(1)

where HSOC is the spin–orbit coupling operator, and the other notation is defined in Table I. The SOC strengths are the quantities that induce population transfer. A similar quantity to Eq. (1) is defined for the S1–T1 spin–orbit coupling, but it is zero by symmetry for planar C2v geometries, and it is small for other geometries. This is consistent with the El-Sayed selection rule68 that nonradiative transitions from an 1* state to a 3ππ* state (which corresponds to the S1 → T2 case) should be orders of magnitude faster than nonradiative transitions from an 1* state to a 3* state (which corresponds to the S1 → T1 case).

TABLE I.

State ordering used in the MCH representation in SHARC. a

StateLabelSMn
S0 
S1 
T11 −1 
T21 −1 
T10 
T20 
T11 +1 
T21 +1 
StateLabelSMn
S0 
S1 
T11 −1 
T21 −1 
T10 
T20 
T11 +1 
T21 +1 
a

Notation: n is the principal quantum number, M is the magnetic spin quantum number, and S is the spin quantum number.

Four sets of semiclassical trajectory calculations were performed, using the SE, CSDM, TSH, and TSH-EDC methods. All dynamics calculations were carried out with eight states in the MCH representation; the eight states are listed in Table I. The initial geometries and velocities were randomly sampled for each of the semiclassical trajectory calculations from the S0 rotationless ground-vibrational-state Wigner distribution69 with MP2 frequencies. The trajectories are propagated using the CASSCF(12,10) potential surfaces and couplings. The initial state is taken as the vertically excited S1 state in the MCH representation. We did not use the time-derivative approximation37 in the CSDM and SE calculations.

If one uses the nonadiabatic couplings (NACs) in the space frame (as are provided by most electronic structure programs), SE, CSDM, and TSH all fail to conserve nuclear angular momentum;70 thus, we modified the usual procedures to conserve angular momentum. In particular, for CSDM and SE, the nuclear equations of motion are propagated using projected NACs to ensure the conservation of angular momentum,70 and for TSH and TSH-EDC trajectories, projected NACs are used as the direction to adjust momentum after the trajectory hops. However, to ensure the better conservation of energy, NACs are used to compute the gradients of the potentials in the diagonal representation21 [the NAC needs to be used because they are in the off-diagonal elements in the MCH representation—see Eq. (13) in a recent review;22 therefore, the gradient in the diagonal representation needs the NACs in the MCH representation]. This step may reintroduce some nonconservation of angular momentum; we examine the question of angular momentum conservation further in the  Appendix.

The nuclear equations of motion were integrated by the velocity-Verlet algorithm with an adaptive time step size (with an energy threshold of 5 × 10−5 eV and a minimum step size of 1 × 10−4 fs). A substep of 200 is used for propagating the electronic coefficients. The calculation of energy gradients of the diagonal states in the current Sharc approach only considers the contribution from the NAC terms [see the second term of Eq. (13) in Ref. 22 and Eq. (48) in Ref. 36]; it ignores the contributions from the Hellmann–Feynman terms because they are hard to calculate and expected to be small in applications such as the present one that have no heavy atoms and no laser field. Thus, the energy need not be precisely converged even in the limit of a very small step size, but in practice, we found that the overall total energy conservation is good enough for the present results to be meaningful. In particular, Fig. S3 in the supplementary material shows the time evolution of the total energy; the total energy is conserved within 80 cm−1 for CSDM, SE, and TSH without decoherence and for most of the TSH-EDC trajectories.

For each dynamics simulation, 200 trajectories were performed with a simulation time of 500 fs for each trajectory.

We consider two singlet states (S0 and S1) and two triplet states (T1 and T2). In the Franck–Condon region, the dominant characters of the S0, S1, T1, and T2 states are 1n2, 1*, 3*, and 3ππ*, respectively. The active orbitals at the ground-state equilibrium structure are shown in Fig. 1. As listed in Table II, the dominant configurations of the SA-CASSCF wave functions for S0, S1, T1, and T2 states at the S0 equilibrium geometry have weights in the 86%–96% range. Since a molecule is usually considered to be multireference when the weight of the dominant configuration is less than or equal to 95%,71 this shows that multireference methods should be used for the present CH2S molecule, which was also demonstrated dynamically by Mai et al.Table II also shows the dominant-configuration weights for the geometry with RC–S = 2.10 Å, and we see that the weights drop to the 62%–93% range, showing even stronger multireference character.

TABLE II.

The characters and the dominant configurations of the states, their dominant configurations, the percentage weight, and the M diagnostics of the dominant states as calculated by CASSCF(12,10) at the S0 equilibrium geometry and at a geometry with a stretched C–S bond.

S0 equilibriumaRC–S = 2.10 Å
PercentageMPercentageM
StateCharacterDominant configurationaweightbdiagnosticcweightbdiagnosticc
S0 1n2 222 222 000 0 85.6 0.10 62.4 0.27 
S1 1* 222 221 100 0 92.4 0.11 87.4 0.13 
T1 3* 222 221 100 0 94.5 0.03 89.9 0.08 
T2 3ππ* 222 212 100 0 96.0 0.03 92.7 0.09 
S0 equilibriumaRC–S = 2.10 Å
PercentageMPercentageM
StateCharacterDominant configurationaweightbdiagnosticcweightbdiagnosticc
S0 1n2 222 222 000 0 85.6 0.10 62.4 0.27 
S1 1* 222 221 100 0 92.4 0.11 87.4 0.13 
T1 3* 222 221 100 0 94.5 0.03 89.9 0.08 
T2 3ππ* 222 212 100 0 96.0 0.03 92.7 0.09 
a

The orbital occupation numbers of the dominant configuration, with the orbitals ordered as in Fig. 1. The occupation numbers of the dominant configurations are the same for the two geometries in this table. The orbitals in active space for RC–S = 2.10 Å (with other geometrical parameters the same as those for the S0 equilibrium structure) are given in Fig. S1 of the supplementary material; they are similar to the ones in Fig. 1.

b

The percentage weight is 100C02, where C0 is the coefficient of the dominant configuration.

c

The M diagnostics are calculated from the state-specific natural-orbital occupation numbers of the state-averaged calculations.

Another way to assess the multiconfigurational character is the M diagnostic, which can be used to label multireference character as small (<0.05), modest (0.05–0.10), or large (>0.10).72Table II shows that the M diagnostic also indicates increasing multireference character as the C–S bond is stretched, and at RC–S = 2.10 Å, two of the states have large multireference character and two have modest multireference character.

The vertical excitation energies calculated with different methods are given in Table III, where they are compared to MS-CASPT2 and MRCISD+P calculations from Ref. 23 (both carried out with the ANO-RCC-VQZP73 basis set) and to the best theoretical estimates of Loos et al.,74 which are estimated to be accurate within ∼0.03 eV. Table III shows that all the CASSCF calculations overestimate the S1 and T1 excitation energies. The CASSCF(12,10) calculation predicts the T2 excitation energy quite accurately compared with the benchmark results, while the CASSCF(10,6) calculation underestimates the T2 excitation energy. The T2–S1 energy gap is known to be critical for the efficiency of ISC in H2CS,23 and this gap is 0.2 eV and 0.4 eV, too small for the large (12,10) and small (12,6) active spaces, respectively.

TABLE III.

Vertical excitation energies and T2–S1 gap (in eV).a

MethodBasis setS1 (1*)T1 (3*)T2 (3ππ*)T2–S1
CASSCF(12,10) 6-31G* 2.43 2.22 3.44 1.01 
CASSCF(10,6) 6-31G* 2.29 2.04 3.08 0.80 
CASSCF(10,6) cc-pVDZ 2.31 2.03 3.09 0.78 
MS-CASPT2(12,10) ANO-RCC-VQZP 2.25 2.00 3.45 1.23 
MRCISD+P(12,10) ANO-RCC-VQZP 2.19 1.91 3.42 1.20 
Theoretical best estimate 2.22 1.94 3.43 1.21 
MethodBasis setS1 (1*)T1 (3*)T2 (3ππ*)T2–S1
CASSCF(12,10) 6-31G* 2.43 2.22 3.44 1.01 
CASSCF(10,6) 6-31G* 2.29 2.04 3.08 0.80 
CASSCF(10,6) cc-pVDZ 2.31 2.03 3.09 0.78 
MS-CASPT2(12,10) ANO-RCC-VQZP 2.25 2.00 3.45 1.23 
MRCISD+P(12,10) ANO-RCC-VQZP 2.19 1.91 3.42 1.20 
Theoretical best estimate 2.22 1.94 3.43 1.21 
a

The results in the first two rows are from the present study, the next three rows are from Mai et al.,23 and the last row is the theoretical best estimate of Loos et al.74 

The dynamical simulation of ISC requires potential energy surfaces not just in the Franck–Condon region. A previous study75 has shown that the C–S bond stretch is the dominant degree of freedom that promotes intersystem crossing in CH2S. The CASSCF(12,10) potential energy curves of the four involved states (S0, S1, T1, and T2 in the MCH representation) are shown in Fig. 2 as functions of the C–S distance. Figure 2 also shows the S1–T2 SOC strength defined in Eq. (1). Figure 2 shows that the energy gaps between the singlet and triplet states are shrunk and the S1–T2 SOC strength increases as the C–S bond is stretched; both of these trends lead to greater ISC at larger RC–S. Eventually, when the potential energy curves have risen to ∼5.0 eV, where the C–S bond is essentially dissociated, the MCH potential curves become nearly degenerate, corresponding to triplet CH2 + triplet S. The S1–T1 crossing occurs at RC–S = 2.05 Å, which is similar to what was calculated with MS-CASPT2 in Ref. 23, as can be seen from Fig. S2 of the supplementary material, especially when one considers that crossing locations are quite sensitive to the electronic structure methods; therefore, the CASSCF calculations with a large active space are quite encouraging for the present CH2S system.

FIG. 2.

Potential energy curves (left scale) and S1–T2 SOC strength (right scale) as functions of the C–S distance from CASSCF(12,10) calculations.

FIG. 2.

Potential energy curves (left scale) and S1–T2 SOC strength (right scale) as functions of the C–S distance from CASSCF(12,10) calculations.

Close modal

Mai et al.23 showed that methods may give state ordering and energies of the low-lying singlet and triplet electronic excited states that are “within an acceptable error of 0.2 eV–0.3 eV,” but still predict qualitatively wrong dynamics. In particular, they found that some dynamical calculations yielded very small triplet populations on the 500 fs time scale, whereas others showed a steady increase in the triplet populations across the time period of the simulations, leading to 3% or more triplet population in that time period. Although the ISC rate of CH2S has apparently not been measured, they concluded that it is small based on the high measured51–53 fluorescence yield and, therefore, that the simulations with a low triplet buildup are the more accurate ones. Next, we look at the triplet yields observed in the present simulations.

The total energy during trajectory propagation in the present study is conserved within 80 cm−1 for CSDM, SE, and TSH without decoherence and most of the TSH-EDC trajectories, as can be seen from Fig. S3 of the supplementary material. Before discussing ensemble-averaged results, we first discuss some trajectories that violate the El-Sayed selection discussed in Sec. II B. We found that 2 out of the 200 trajectories in the CSDM simulations showed population transfer from the S1 state to the T1 state followed by decohering (pointer-state switching) to the T1 state, as shown in the left panel of Fig. S4 of the supplementary material. Further analysis shows that the switch of the pointer state from S1 to T1 in the two CSDM trajectories happens when the switching probability (to T1) is on the order of 1 × 10−6 with an even smaller generated random number, as indicated in Table S1 of the supplementary material. 13 out of 200 trajectories in the TSH-EDC simulations have a similar population transfer (three typical trajectories are depicted in the right panel of Fig. S4). Of course, the El-Sayed selection rule, like most selection rules, is really just a propensity rule; trajectories that violate it are not actually forbidden; they are just expected to be rare. Since our goal here is to analyze the major pathways and not have the conclusions muddied by contributions from rare events (because it would require larger ensembles to include such rare events in a statistically meaningful way), we eliminated these anti-El-Sayed trajectories from the ensemble averages. Therefore, 198, 200, 187, and 200 trajectories are included in the ensemble-average analysis for CSDM, SE, TSH-EDC, and TSH simulations, respectively. We will come back to these rare events later.

The ensemble-average time evolution of the state populations for the two triplet states and the SOC strengths between S1 and two triplet states are shown in Fig. 3 for the four ensembles of trajectories (CSDM, SE, TSH-EDC, and TSH). The state populations given in Fig. 3 are calculated from the squares of the state coefficients in the spin-free MCH representation.22 We first analyze the simulations carried out by the two methods that include decoherence, namely, the CSDM and TSH-EDC trajectories that are depicted in Figs. 3(a) and 3(c). These figures show that similar results are obtained for these two quite different types of nonadiabatic trajectory methods. The time evaluation of the T2 population, which is shown as a solid orange curve, shows a damped oscillation. The oscillation period perfectly matches the variation of the spin–orbit couplings, which vary due to the oscillation of the C–S distance. The ensemble-averaged C–S distance as a function of time is shown in Fig. 4.

FIG. 3.

Ensemble-averaged time evolution of the state populations for T1 and T2 excited states in MCH representation and the SOC strengths for S1–T1 and S1–T2. (a) CSDM, (b) SE, (c) TSH-EDC, and (d) TSH without decoherence. The state populations are the solid curves with the scale on the left; the SOC strengths are the dotted curves with the scales on the right; the right-hand scale has a break because the SOC strength is much larger for S1–T2 than that for S1–T1. The numbers on the ordinate scale must be multiplied by the factors in parentheses in the ordinate label; thus, for example, the tick mark at 0.2 denotes 0.2 × 10−2, i.e., 0.2%.

FIG. 3.

Ensemble-averaged time evolution of the state populations for T1 and T2 excited states in MCH representation and the SOC strengths for S1–T1 and S1–T2. (a) CSDM, (b) SE, (c) TSH-EDC, and (d) TSH without decoherence. The state populations are the solid curves with the scale on the left; the SOC strengths are the dotted curves with the scales on the right; the right-hand scale has a break because the SOC strength is much larger for S1–T2 than that for S1–T1. The numbers on the ordinate scale must be multiplied by the factors in parentheses in the ordinate label; thus, for example, the tick mark at 0.2 denotes 0.2 × 10−2, i.e., 0.2%.

Close modal
FIG. 4.

The evolution of the C–S distance with (a) CSDM, (b) SE, (c) TSH-EDC, and (d) TSH without decoherence. Individual results are plotted in pink, and the ensemble-averaged result is in blue.

FIG. 4.

The evolution of the C–S distance with (a) CSDM, (b) SE, (c) TSH-EDC, and (d) TSH without decoherence. Individual results are plotted in pink, and the ensemble-averaged result is in blue.

Close modal

The damped oscillation curves for CSDM and TSH-EDC indicate that the T2 population tends to about 0.2% (after a sufficiently long simulation time). A similar oscillation in the T2 population was obtained by Mai et al.23 in TSH-EDC dynamics based on MS-CASPT2(10,6)23 electronic structure calculations, although their reported triplet population was smaller, oscillating between 0% and 0.2% (as read from their Fig. 2); based on this observation, Mai et al.23 concluded that no persistent triplet population builds up on a subpicosecond simulation time scale. Our trajectories do show some pointer-state switching (in CSDM) and surface hopping (in TSH-EDC) from S1 to triplet states, but the trend—such as that in Ref. 23—is that the amount of triplet state population is small and steady [as compared to the increasing populations in Figs. 3(b) and 3(d), which are discussed below].

Mai et al.23 also inferred that the oscillations are mainly due to use of the spin-mixed potential energy surfaces (the diagonal representation) for the propagation of intersystem crossing trajectories in Sharc. The oscillation period extracted from the variation of the T2 population, of the S1–T2 SOC strength, and of the C–S distance in Figs. 3(a), 3(c), 4(a), and 4(c) is about 46 fs (with about eleven cycles in 500 fs simulation time). This is in good agreement with simulation results23 carried out by MS-CASPT2, which showed an oscillation period of 43 fs. These periods correspond to frequencies of 725 cm−1 and 775 cm−1, respectively, which are in good agreement with the observed frequency51 for C–S stretching in the S1 state of CH2S, which is 820 cm−1.

The current CASSCF(12,10) results for the triplet populations differ from another set of TSH-EDC simulations in Ref. 23 that used SA-CASSCF(10,6) and from the generalized ab initio multiple spawning (GAIMS) simulations of Ref. 75 with an even smaller (4,3) active space; both the results with a small active space lead to a constant increase in the T2 population. Figure S2 of the supplementary material compares potential energy curves obtained with different electronic structure methods, and it shows that the CASSCF results with a large active space in the present work are more in line with the MS-CASPT2(10,6) ones. This clearly shows the critical role of the active space in photodynamics simulations. It is encouraging that the SA-CASSCF method with the larger active space used here can lead to reasonable results (the inference that it is reasonable is based on comparing to the results with MS-CASPT2,23 which are consistent with experimental findings51–53), whereas the small-active-space CASSCF calculations of Ref. 23 does not.

We draw two main conclusions from the CSDM and TSH-EDC calculations. (1) With a well-chosen active space, dynamical results calculated with CASSCF(10,12) agree well with previous results using MS-CASPT2, which is a more expensive method whose cost is sometimes unaffordable. (2) It is very encouraging that the simple TSH-EDC calculations agree so well with the CSDM calculations, which have been shown in past work to be more accurate than TSH calculations without the EDC corrections.43,76,77

Next, we consider the simulations carried out without decoherence, in particular, the SE and TSH results in Figs. 3(b) and 3(d). These two methods give very similar results to one another for the present tests (which is not always the case for other systems), but the results are strikingly different from CSDM and TSH-EDC. Whereas CSDM shows T2 oscillating around 0.2% and T1 oscillating around 0.003%, the calculations without decoherence show T2 oscillating around 0.4% and T1 gradually increasing to 0.5%, and still increasing at the end of the simulation at 500 fs. That is, the T1 populations are not damped, which is qualitatively different from the results including decoherence. Figure 3 shows that the ensemble averages of the spin–orbit coupling strengths are similar for the calculations with and without decoherence, indicating that the trajectories reach similar regions of space. Thus, the differences are due, not so much, to different nuclear trajectories as to the way that decoherence is treated. This is an illustration of the importance of including decoherence in nonadiabatic trajectory simulations. (A similar effect of decoherence was seen in Ref. 47.)

The results obtained here without decoherence are qualitatively similar to earlier dynamics calculations75 with generalized ab initio multiple spawning (GAIMS), in which the triplet population also increased without damping, in particular, to 11% by 150 fs and to 14% by 200 fs. This might come from the smaller active space used in GAIMS, and it would be interesting to do more direct comparisons.

Figure 5 shows ensemble averaged energies of the states; it shows damped oscillations that correlate with the damped oscillations in the C–S distances that are shown in Fig. 4.

FIG. 5.

Time evolution of ensemble-averaged energies of each state (S0, S1, T1, and T2 in the MCH representation) during the 500 fs dynamics simulation with (a) CSDM, (b) SE, (c) TSH-EDC, and (d) TSH without decoherence. In (a) and (b), the plots also show the ensemble average of the self-consistent potential.

FIG. 5.

Time evolution of ensemble-averaged energies of each state (S0, S1, T1, and T2 in the MCH representation) during the 500 fs dynamics simulation with (a) CSDM, (b) SE, (c) TSH-EDC, and (d) TSH without decoherence. In (a) and (b), the plots also show the ensemble average of the self-consistent potential.

Close modal

Figure S5 of the supplementary material shows the time evolution when we include all 200 trajectories (even the anti-El-Sayed ones) in the ensemble average. Due to the rare events, i.e., switching of the pointer state or hopping of the active state, from S1 to T1, both CSDM and TSH-EDC results exhibit much larger anti-El-Sayed T1 populations than SE or TSH, but the populations are still small in an absolute sense, less than 5%. However, it is a well-known deficiency of Monte Carlo sampling methods that one cannot trust rare-event probabilities quantitatively unless one uses large ensembles or rare event sampling.

In this work, we first calculated the vertical excitation energies and the characters of the two lowest singlets and two lowest triplets of the CH2S molecule by the SA-CASSCF method with a reasonably large active space, 12 electrons in 10 orbitals. We found that calculations with CASSCF(12,10) have good accuracy for both the potential energy surfaces and the singlet–triplet gap.

We then performed intersystem crossing simulations by direct dynamics with CASSCF(12,10) to study the effect of decoherence on the dynamics. Four nonadiabatic dynamics methods, namely, CSDM, SE, TSH with an EDC correction, and TSH without a decoherence correction, have been tested. We found that the single-surface propagation of TSH with energy-based decoherence and the self-consistent potential propagation of coherent switching with decay of mixing give virtually identical results even though they incorporate decoherence in quite different ways. The state populations yielded by the two methods with decoherence included agree qualitatively with the state populations obtained by Mai et al.23 with TSH-EDC, even though our calculations are based on CASSCF and their calculations are based on more expensive MS-CASPT2 potentials. The good agreement that we observe between CSDM and TSH-EDC is especially striking when compared to the qualitatively different behavior we see when decoherence is not included, which is the case in conventional surface hopping and the semiclassical Ehrenfest method. Whereas the triplet state populations show damped oscillations toward a small limiting value when decoherence is included, they show steadily increasing behavior when density matrix decoherence is not included. This comparison illustrates the importance of taking electronic decoherence into account in semiclassical nonadiabatic dynamics.

See the supplementary material for active orbitals at the geometry with a stretched C–H bond, representative trajectories with unusual state-population histories, the time evolution of average state populations when these are included in the ensemble averages, coordinates and harmonic frequencies of the ground-state equilibrium geometry, and the template control file used in the CSDM calculations with SHARC.

L.Z. and Y.S. contributed equally to this work.

This work was supported, in part, by the National Natural Science Foundation of China (Grant No. 51536002) and by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences (Award No. DE-SC0015997).

The authors declare no competing financial interest.

The data that support the findings of this study are available within this article and its supplementary material.

In a previous study,70 we have shown that the use of space-frame NACs in direct nonadiabatic dynamics algorithms will not conserve the total angular momentum because in electronic structure codes, the NACs are calculated by varying nuclear Cartesian coordinates. We also showed that the use of projected NACs can eliminate the terms causing this nonconservation problem. In treating ISC, a new term (HSOC) is added to the Hamiltonian, and the issue of nonconservation of nuclear angular momentum needs be revisited because the correction of the gradients of the states in a diagonal representation by the space-frame NACs introduces a new mechanism to violate conservation of nuclear angular momentum.

Figure 6 shows the average time evolution of the angular momentum magnitude |J| along the trajectories for all four test sets. It can be seen that with the NAC corrections of gradients in the diagonal representation, the angular momentum is not strictly conserved as the trajectory propagates. The intersystem crossing trajectories propagated on self-consistent potentials show slightly better conservation than the ones propagated on single-state potentials. The overall nonconvergence caused by the additional HSOC term is not as severe as seen in our previous tests when unprojected NACs were used, and it is not a serious problem, especially when one considers that changes of angular momentum less than 12 are not physically meaningful in semiclassical calculations.

FIG. 6.

Time evolution of the angular momentum magnitude in intersystem crossing dynamics with (a) CSDM, (b) SE, (c) TSH-EDC, and (d) TSH without decoherence. Individual trajectories are shown in pink, and the ensemble is shown in blue.

FIG. 6.

Time evolution of the angular momentum magnitude in intersystem crossing dynamics with (a) CSDM, (b) SE, (c) TSH-EDC, and (d) TSH without decoherence. Individual trajectories are shown in pink, and the ensemble is shown in blue.

Close modal
1.
J. M.
Simmie
,
Prog. Energy Combust. Sci.
29
,
599
(
2003
).
2.
P.
Zhao
,
W.
Yuan
,
H.
Sun
,
Y.
Li
,
A. P.
Kelley
,
X.
Zheng
, and
C. K.
Law
,
Proc. Combust. Inst.
35
,
309
(
2015
).
3.
W. K.
Metcalfe
,
S.
Dooley
, and
F. L.
Dryer
,
Energy Fuels
25
,
4915
(
2011
).
4.
A.
Caracciolo
,
G.
Vanuzzo
,
N.
Balucani
,
D.
Stranges
,
P.
Casavecchia
,
L.
Pratali Maffei
, and
C.
Cavallotti
,
J. Phys. Chem. A
123
,
9934
(
2019
).
5.
C.
Adachi
,
M. A.
Baldo
,
M. E.
Thompson
, and
S. R.
Forrest
,
J. Appl. Phys.
90
,
5048
(
2001
).
6.
X.
Yang
,
D.
Neher
,
D.
Hertel
, and
T. K.
Däubler
,
Adv. Mater.
16
,
161
(
2004
).
7.
J. L.
Segura
,
N.
Martín
, and
D. M.
Guldi
,
Chem. Soc. Rev.
34
,
31
(
2005
).
8.
J. C.
Deaton
,
S. C.
Switalski
,
D. Y.
Kondakov
,
R. H.
Young
,
T. D.
Pawlik
,
D. J.
Giesen
,
S. B.
Harkins
,
A. J. M.
Miller
,
S. F.
Mickenberg
, and
J. C.
Peters
,
J. Am. Chem. Soc.
132
,
9499
(
2010
).
9.
H.
Yersin
,
A. F.
Rausch
,
R.
Czerwieniec
,
T.
Hofbeck
, and
T.
Fischer
,
Coord. Chem. Rev.
255
,
2622
(
2011
).
10.
H.
Uoyama
,
K.
Goushi
,
K.
Shizu
,
H.
Nomura
, and
C.
Adachi
,
Nature
492
,
234
(
2012
).
11.
T.
Kinoshita
,
J. T.
Dy
,
S.
Uchida
,
T.
Kubo
, and
H.
Segawa
,
Nat. Photonics
7
,
535
(
2013
).
12.
Y.
Shu
and
B. G.
Levine
,
J. Chem. Phys.
142
,
104104
(
2015
).
13.
W.
Chang
,
D. N.
Congreve
,
E.
Hontz
,
M. E.
Bahlke
,
D. P.
McMahon
,
S.
Reineke
,
T. C.
Wu
,
V.
Bulović
,
T.
Van Voorhis
, and
M. A.
Baldo
,
Nat. Commun.
6
,
6415
(
2015
).
14.
M.
Richter
,
P.
Marquetand
,
J.
González-Vázquez
,
I.
Sola
, and
L.
González
,
J. Chem. Phys. Lett.
3
,
3090
(
2012
).
15.
R. A.
Vogt
,
C.
Reichardt
, and
C. E.
Crespo-Hernández
,
J. Phys. Chem. A
117
,
6580
(
2013
).
16.
R. R.
Zaari
and
S. A.
Varganov
,
J. Phys. Chem. A
119
,
1332
(
2015
).
17.
H. A.
Ernst
,
T. J. A.
Wolf
,
O.
Schalk
,
N.
González-García
,
A. E.
Boguslavskiy
,
A.
Stolow
,
M.
Olzmann
, and
A.-N.
Unterreiner
,
J. Phys. Chem. A
119
,
9225
(
2015
).
18.
T.
Schnappinger
,
M.
Marazzi
,
S.
Mai
,
A.
Monari
,
L.
González
, and
R.
de Vivie-Riedle
,
J. Chem. Theory Comput.
14
,
4530
(
2018
).
19.
J. P.
Zobel
and
L.
González
,
ChemPhotoChem
3
,
833
(
2019
).
20.
M.
Richter
,
P.
Marquetand
,
J.
González-Vázquez
,
I.
Sola
, and
L.
González
,
J. Chem. Theory Comput.
7
,
1253
(
2011
).
21.
S.
Mai
,
P.
Marquetand
, and
L.
González
,
Int. J. Quantum Chem.
115
,
1215
(
2015
).
22.
S.
Mai
,
P.
Marquetand
, and
L.
González
,
Wiley Interdiscip. Rev.: Comput. Mol. Sci.
8
,
e1370
(
2018
).
23.
S.
Mai
,
A. J.
Atkins
,
F.
Plasser
, and
L.
González
,
J. Chem. Theory Comput.
15
,
3470
(
2019
).
24.
R.
Valero
and
D. G.
Truhlar
,
J. Phys. Chem. A
111
,
8536
(
2007
).
25.
A.
Bjerre
and
E. E.
Nikitin
,
Chem. Phys. Lett.
1
,
179
(
1967
).
26.
J. C.
Tully
and
R. K.
Preston
,
J. Chem. Phys.
55
,
562
(
1971
).
27.
P. J.
Kuntz
,
J.
Kendrick
, and
W. N.
Whitton
,
Chem. Phys.
38
,
147
(
1979
).
28.
D. G.
Truhlar
,
J. W.
Duff
,
N. C.
Blais
,
J. C.
Tully
, and
B. C.
Garrett
,
J. Chem. Phys.
77
,
764
(
1982
).
29.
S.
Chapman
, in , edited by
M.
Baer
and
C. Y.
Ng
(
John Wiley & Sons, Inc, Hoboken, NJ
,
1992
), Vol. 82, pp.
423
483
.
30.
J. C.
Tully
,
J. Chem. Phys.
93
,
1061
(
1990
).
31.
C.
Zhu
,
S.
Nangia
,
A. W.
Jasper
, and
D. G.
Truhlar
,
J. Chem. Phys.
121
,
7658
(
2004
).
32.
H. D.
Meyer
and
W. H.
Miller
,
J. Chem. Phys.
70
,
3214
(
1979
).
33.
D. A.
Micha
,
J. Chem. Phys.
78
,
7138
(
1983
).
34.
M. S.
Topaler
,
T. C.
Allison
,
D. W.
Schwenke
, and
D. G.
Truhlar
,
J. Chem. Phys.
109
,
3321
(
1998
).
35.
D. A.
Micha
,
J. Phys. Chem. A
103
,
7562
(
1999
).
36.
Y.
Shu
,
L.
Zhang
,
S.
Mai
,
S.
Sun
,
L.
González
, and
D. G.
Truhlar
,
J. Chem. Theory Comput.
16
,
3464
(
2020
).
37.
Y.
Shu
,
L.
Zhang
,
S.
Sun
, and
D. G.
Truhlar
,
J. Chem. Theory Comput.
16
,
4098
(
2020
).
38.
B. J.
Schwartz
and
P. J.
Rossky
,
J. Chem. Phys.
101
,
6902
(
1994
).
39.
E. R.
Bittner
and
P. J.
Rossky
,
J. Chem. Phys.
103
,
8130
(
1995
).
40.
B. J.
Schwartz
,
E. R.
Bittner
,
O. V.
Prezhdo
, and
P. J.
Rossky
,
J. Chem. Phys.
104
,
5942
(
1996
).
41.
M. D.
Hack
and
D. G.
Truhlar
,
J. Chem. Phys.
114
,
9305
(
2001
).
42.
C.
Zhu
,
A. W.
Jasper
, and
D. G.
Truhlar
,
J. Chem. Phys.
120
,
5543
(
2004
).
43.
C.
Zhu
,
A. W.
Jasper
, and
D. G.
Truhlar
,
J. Chem. Theory Comput.
1
,
527
(
2005
).
44.
G.
Granucci
and
M.
Persico
,
Theor. Chem. Acc.
117
,
1131
(
2007
).
45.
G.
Granucci
,
M.
Persico
, and
A.
Zoccante
,
J. Chem. Phys.
133
,
134111
(
2010
).
46.
F.
Plasser
,
S.
Mai
,
M.
Fumanal
,
E.
Gindensperger
,
C.
Daniel
, and
L.
González
,
J. Chem. Theory Comput.
15
,
5031
(
2017
).
47.
A. W.
Jasper
and
D. G.
Truhlar
,
J. Chem. Phys.
127
,
194306
(
2007
).
48.
B. O.
Roos
,
P. R.
Taylor
, and
P. E. M.
Sigbahn
,
Chem. Phys.
48
,
157
(
1980
).
49.
J.
Finley
,
P.-Å.
Malmqvist
,
B. O.
Roos
, and
L.
Serrano-Andrés
,
Chem. Phys. Lett.
288
,
299
(
1998
).
50.
P.
Celani
and
H.-J.
Werner
,
J. Chem. Phys.
112
,
5546
(
2000
).
51.
D. J.
Clouthier
and
D. A.
Ramsay
,
Annu. Rev. Phys. Chem.
34
,
31
(
1983
).
52.
M.
Kawasaki
,
K.
Kasatani
,
Y.
Ogawa
, and
H.
Sato
,
Chem. Phys.
74
,
83
(
1983
).
53.
D. C.
Moule
and
E. C.
Lim
,
J. Phys. Chem. A
106
,
3072
(
2002
).
54.
W. J.
Hehre
,
L.
Radom
,
P. V. R.
Schleyer
, and
J. A.
Pople
,
Ab Initio Molecular Orbital Theory
(
Wiley
,
New York
,
1986
).
55.
H.-J.
Werner
,
P. J.
Knowles
,
G.
Knizia
,
F. R.
Manby
, and
M.
Schütz
,
Wiley Interdiscip. Rev.: Comput. Mol. Sci.
2
,
242
(
2012
).
56.
H.-J.
Werner
,
P. J.
Knowles
,
F. R.
Manby
,
J. A.
Black
,
K.
Doll
,
A.
Heßelmann
,
D.
Kats
,
A.
Köhn
,
T.
Korona
,
D. A.
Kreplin
,
Q.
Ma
,
T. F.
Miller
,
A.
Mitrushchenkov
,
K. A.
Peterson
,
I.
Polyak
,
G.
Rauhut
, and
M.
Sibaev
,
J. Chem. Phys.
152
,
144107
(
2020
).
57.
F.
Aquilante
,
J.
Autschbach
,
R. K.
Carlson
,
L. F.
Chibotaru
,
M. G.
Delcey
,
L.
De Vico
,
I. F.
dez Galván
,
N.
Ferré
,
L. M.
Frutos
,
L.
Gagliardi
,
M.
Garavelli
,
A.
Giussani
,
C. E.
Hoyer
,
G.
Li Manni
,
H.
Lischka
,
D.
Ma
,
P. Å.
Malmqvist
,
T.
Müller
,
A.
Nenov
,
M.
Olivucci
,
T. B.
Pedersen
,
D.
Peng
,
F.
Plasser
,
B.
Pritchard
,
M.
Reiher
,
I.
Rivalta
,
I.
Schapiro
,
J.
Segarra-Martí
,
M.
Stenrup
,
D. G.
Truhlar
,
L.
Ungur
,
A.
Valentini
,
S.
Vancoillie
,
V.
Veryazov
,
V. P.
Vysotskiy
,
O.
Weingart
,
F.
Zapata
, and
R.
Lindh
,
J. Comput. Chem.
37
,
506
(
2016
).
58.
T.
Lu
and
F.
Chen
,
J. Comput. Chem.
33
,
580
(
2012
).
59.
W.
Humphrey
,
A.
Dalke
, and
K.
Schulten
,
J. Mol. Graphics
14
,
33
(
1996
).
60.
J. A.
Karwowski
and
I.
Shavitt
, in
Handbook of Molecular Physics and Quantum Chemistry
, edited by
S.
Wilson
(
John Wiley & Sons
,
Chichester
,
2003
), p.
227
.
61.
J. A.
Pople
,
R.
Seeger
, and
R.
Krishnan
,
Int. J. Quantum Chem.
12
,
149
(
1977
).
62.
Y.
Shu
,
L.
Zhang
, and
D. G.
Truhlar
, SHARC-MN, version 1.0,
University of Minnesota
,
Minneapolis
,
2020
, https://comp.chem.umn.edu/sharcmn.
63.
S.
Mai
,
M.
Richter
,
M.
Heindl
,
M. F. S. J.
Menger
,
A.
Atkins
,
M.
Ruckenbauer
,
F.
Plasser
,
L. M.
Ibele
,
S.
Kropf
,
M.
Oppel
,
P.
Marquetand
, and
L.
González
, SHARC2.1,
University of Vienna
,
Wien
,
2019
, https://sharc-md.org.
64.
P. C.
Hariharan
and
J. A.
Pople
,
Theor. Chim. Acta
28
,
213
(
1973
).
65.
V. A.
Rassolov
,
M. A.
Ratner
,
J. A.
Pople
,
P. C.
Redfern
, and
L. A.
Curtiss
,
J. Comput. Chem.
22
,
976
(
2001
).
66.
M.
Reiher
and
A.
Wolf
,
J. Chem. Phys.
121
,
2037
(
2004
).
67.
M.
Reiher
and
A.
Wolf
,
J. Chem. Phys.
121
,
10945
(
2004
).
68.
M. A.
El-Sayed
,
Acc. Chem. Res.
1
,
8
(
1968
).
70.
Y.
Shu
,
L.
Zhang
,
Z.
Varga
,
K. A.
Parker
,
S.
Kanchanakungwankul
,
S.
Sun
, and
D. G.
Truhlar
,
J. Phys. Chem. Lett.
11
,
1135
(
2020
).
71.
S. O.
Odoh
,
G. L.
Manni
,
R. K.
Carlson
,
D. G.
Truhlar
, and
L.
Gagliardi
,
Chem. Sci.
7
,
2399
(
2016
).
72.
O.
Tishchenko
,
J.
Zheng
, and
D. G.
Truhlar
,
J. Chem. Theory Comput.
4
,
1208
(
2008
).
73.
B. O.
Roos
,
R.
Lindh
,
P.-Å.
Malmqvist
,
V.
Veryazov
, and
P.-O.
Widmark
,
J. Phys. Chem. A
109
,
6575
(
2005
).
74.
P.-F.
Loos
,
A.
Scemama
,
A.
Blondel
,
Y.
Garniron
,
M.
Caffarel
, and
D.
Jacquemin
,
J. Chem. Theory Comput.
14
,
4360
(
2018
).
75.
B. F. E.
Curchod
,
C.
Rauer
,
P.
Marquetand
,
L.
González
, and
T. J.
Martínez
,
J. Chem. Phys.
144
,
101102
(
2016
).
76.
A. W.
Jasper
,
S.
Nangia
,
C.
Zhu
, and
D. G.
Truhlar
,
Acc. Chem. Res.
39
,
101
(
2006
).
77.
D. G.
Truhlar
, in
Quantum Dynamics of Complex Molecular Systems
, edited by
D. A.
Micha
and
I.
Burghardt
(
Springer
,
Berlin
,
2007
), pp.
227
243
.

Supplementary Material