Plain language summary
When a molecule absorbs light, its electrons are excited. The molecular geometry distorts until the excess of energy dissipates as heat or is reemitted as light. It is possible to simulate using a computer how a molecule reacts to light excitation. Still, such dynamics simulations involve complex and slow quantum-mechanical methods. In this work, we present a new approach to compute one of the quantities needed for the simulations, the nonadiabatic couplings. Based on a new theory proposed by Baeck and An, our method (we call it TD-BA) gets the couplings from electronic energies and their variations with time. Because these quantities are always available during dynamics simulations, TD-BA couplings are much simpler and faster to compute than usual couplings. We tested our TD-BA method for the simulation of two molecules, ethylene and fulvene. The results show that TD-BA does a good job predicting the dynamics of these molecules.
I. Introduction
Surface hopping has become one of the main tools for nonadiabatic dynamics simulations of photoexcited molecules1,2. Its popularity stems from its underlying independent trajectory and local approximations, requiring electronic quantities computed only at the nuclear geometry of classical trajectories. Moreover, surface hopping is usually done in adiabatic representation, avoiding cumbersome diabatization procedures usually required for wavepacket propagation. Thus, surface hopping is commonly implemented to work with on-the-fly calculations of electronic properties (energies, energy gradients, and couplings) obtained with many different quantum chemical methods3–9, and it requires neither pre-calculated nor modeled global potential energy surfaces.
The primary constraint for interfacing surface hopping to a new quantum chemical method or program is to get the first-order nonadiabatic coupling vectors between adiabatic electronic states J and L. These quantities are not as widespread in quantum-chemical packages as excited-state energy gradients for instance, though the situation has improved in the last two decades10–12. Alternatively, the coupling vectors can also be obtained from an approximation involving potential energy Hessian matrices13,14. Moreover, there are surface hopping variants that do not even require explicit nonadiabatic coupling calculations, as is the case of the Zhu-Nakamura15,16 and Belyaev-Lebedev17–19 surface hopping approaches. The former estimates nonadiabatic events from minimum energy gaps and energy gradients around this minimum, while the latter does from minimum energy gaps and second time-derivative of the gap at the minimum gap. The popular fewest switches surface hopping (FSSH)20 requires nonadiabatic coupling vectors. Nonetheless, because these vectors are always projected on the velocities within the formalism, they can be replaced by time-derivative couplings , which are much simpler to implement than nonadiabatic coupling vectors. In general, the calculation of time-derivative couplings follows the Hammes-Schiffer-Tully prescription21, and they are obtained from wavefunction overlaps22–24. Alternatively, the coupling vector calculation can be avoided entirely in FSSH, by employing the local diabatization formalism25, which is also based on wavefunction overlaps.
All these models possess algorithmic limitations for their universal adoption. The Hessian-based nonadiabatic coupling approximation is too costly for on-the-fly dynamics, although it has been used with machine-learning (ML) nonadiabatic dynamics26. Zhu-Nakamura and Belyaev-Lebedev approaches require determination of the diabatic crossing point, which, translated into adiabatic representation, means that the hopping at time t depends not only on the previous timesteps but also on the future of the trajectory. (These methods are usually implemented in terms of a three-point approximation15,17–19, evaluating the hopping at t from information computed at t – Δt, t, and t + Δt.) On the other hand, the wavefunction overlaps required by time-derivative coupling and local diabatization approaches boil down to nonorthogonal atomic-orbital overlap matrices, demanding lengthy algorithmic intervention whenever a new method or program needs to be interfaced to surface hopping.
Recently, Baeck and An27 conjectured that nonadiabatic couplings near an energy crossing point could be estimated from the simple 1D model
involving only the energy gap ΔEJL = EJ – EL and the second derivative of this gap with respect to the coupling direction Q. Details of their model are surveyed in Section II. Baeck and An’s numerical tests for different molecules have indeed shown an outstanding agreement between this simple approximation and the exact nonadiabatic coupling27,28. Our goal in this paper is to expand and explore the Baeck-An (BA) model to be used in surface hopping dynamics.
Baeck and An pointed out two potential limitations of their model27. First, it is strictly valid for molecular configurations coupling only two states. Second, the model is restricted to a single degree of freedom. When using their model for dynamics, there is nothing that we can do about the first limitation but to hope that two-state coupling dominates the nonadiabatic interactions. However, concerning the second limitation, it is fortunate that FSSH does not require the full nonadiabatic coupling vector but only its projection on the velocity direction. As we shall discuss, we can explore this feature to expand the 1D BA model to multidimensional dynamics, using a time-dependent (TD) BA model.
The TD-BA model is not supposed to replace the exact calculation of nonadiabatic couplings in surface hopping. It can be, however, an alternative approach to be used with electronic structure methods for which nonadiabatic couplings are still unavailable or are too computationally expensive to be computed on-the-fly. The TD-BA model can also be used as an inexpensive diagnostic of nonadiabatic interactions during exploratory dynamics, triggering higher-level calculations. Another potential use of the TD-BA model is in ML. Several authors have reported difficulties in the ML-prediction of nonadiabatic couplings29–31. Their narrow functional shapes require too large training sets and active learning, being the main challenge for developing ML nonadiabatic dynamics. The TD-BA model, requiring only predicting energy gaps, may open a new avenue to overcome these problems.
In this work we show how TD-BA couplings can be seamlessly integrated into FSSH. We tested it by performing decoherence-corrected (DC) FSSH dynamics for ethylene and fulvene (Figure 1) using TD-BA couplings and comparing the results to DC-FSSH dynamics of these molecules using full nonadiabatic coupling vectors. Dynamics of both systems have been extensively studied before, and, recently, Ibele and Curchod32 proposed that they may be considered molecular examples for two of the three basic 1D Tully models20. Ethylene is a crucial prototype32 for internal conversion in various systems as diverse as protonated Schiff bases33, substituted polyenes34, and purine bases35. In turn, fulvene poses the challenge of S1↔S0 recurrences at an extended crossing seam with both peaked and sloped conical intersections36. The results for both molecules are encouraging, showing that TD-BA-based dynamics delivers a semiquantitative agreement with reference datasets computed with exact couplings.
Figure 1. Molecular structure of ethylene and fulvene.
II. The BA model for nonadiabatic couplings
Baeck and An have shown that the nonadiabatic coupling maximum for a two-states system can be approximated as27
where Q is a one-dimensional coordinate coupling the adiabatic states J and L. Qc is the geometry where the magnitude of the coupling is the largest. This expression arises from assuming that the nonadiabatic coupling can be represented as a Lorenz function of Q,
where κ, the coupling constant, and Δc, the adiabatic energy gap at Qc, are constants. At Qc, the maximum nonadiabatic coupling is simply
Δc, which by construction corresponds to twice the diabatic coupling, is given as the adiabatic energy gap between J and L at Qc
In turn, to get the coupling constant κ, Baeck and An explore the similarities between their model and the linear vibronic coupling (LVC) model37. While in LVC, Köppel, Gronki, and Mahapatra propose to compute κ from
Baeck and An propose the expression
which, compared to LVC, equivales to introduce an asymmetric linear term in Q in the expansion of ΔEJL. This term aims to make the model more flexible so that the minimum of ΔEJL does not necessarily coincide with Qc.
Finally, the nonadiabatic coupling maximum in Equation (2) is obtained by replacing κ given in Equation (7) into Equation (4).
Baeck and An also conjectured that the functional form of Equation (2) extended to an arbitrary point Q should describe the nonadiabatic coupling also in the region near Qc. This is the source of Equation (1).
III. TD-BA couplings
To employ the BA model in surface hopping, we start by writing their expression as a vector in the unitary direction of Q
where R is the full-dimensional molecular geometry, which has a crossing at Rc. In the spirit of the Baeck-An 1D-approximation, Equation (8) assumes that the only non-null energy Hessian component is that one in the Q direction.
Using the chain rule, we can rewrite the Hessian as a time derivative
where vQ is the projection of the nuclear velocity on . Replacing Equation (9) in Equation (8) yields the time-dependent version of the Baeck-An nonadiabatic coupling
which implicitly assumes that direction is constant. In the previous expression, tc is the time at which the molecule reaches Rc.
In FSSH20, the quantity of interest is the time-derivative coupling (see Section IV), which using Equation (10) simplifies to
Note that all references to Q cancel out in Equation (11).
To enforce the condition t ≈ tc, we notice that whenever the molecule moves near a minimum of ΔEJL, the square-root argument in Equation (11) must be positive. Thus, we assume that the coupling is null for negative arguments, and the TD-BA nonadiabatic coupling becomes
Although the conditions in Equation (12) are not sufficient to detect t ≈ tc (the square root argument may still be positive away from the minimum gap), they make a good job screening the relevant crossing regions during dynamics, as our tests show.
IV. Using the TD-BA model in surface hopping
The FSSH hopping probability from L to J is computed as20
where σLJ is the time derivative coupling given by the TD-BA approximation in Equation (12), and the complex coefficients c are the solution of the locally-approximated time-dependent electronic Schrödinger equation
Δτ in Equation (13) is the timestep for integrating Equation (14). It is, in general, much smaller than Δt, the timestep to integrate the classical equations of motion.
As usual, the FSSH over-coherence should be corrected and, here, we adopted the simplified decay of mixing38.
Because the h direction is unknown, the velocity adjustment after hopping should be made in the direction of the linear momentum p39.
V. Computational details and datasets
We have used ethylene and fulvene (Figure 1) dynamics to test the method. All calculations were done with the multiconfigurational self-consistent field (MCSCF) method state-averaged over three states for ethylene and two states for fulvene. Details on the computation level of ethylene are discussed in Ref. 39 In brief, the active space was composed of five separate complete active subspaces: the first subspace contained four electrons and four orbitals (π, σ, π*, σ*), and each of the other four subspaces contained two electrons in two orbitals (σ, σ*). Thus, the active space can be denoted as [CAS(4,4) ⊕ 4×CAS(2,2)]. For fulvene, we adopted a standard CAS(6,6) complete active space. The 6-31G(d) basis set40 was used in both cases.
Dynamics was simulated with decoherence-corrected38 fewest switches surface hopping20 for a maximum of 200 fs for ethylene and 60 fs for fulvene. The classical equation were integrated with Δt = 0.1 fs steps. The quantum equations were integrated with Δτ = 0.005 fs steps, using interpolated electronic quantities between classical steps. Decoherence corrections were applied with the standard 0.1 au parameter38.
For ethylene, one set of 500 trajectories was run with the TD-BA model. For the analysis, the TD-BA dataset was compared to surface hopping based on exact nonadiabatic coupling vectors with h-adjusted and p-adjusted velocities. These two datasets, also composed of 500 trajectories each, are the same ones discussed in Ref. 39. For fulvene, we run three datasets of 200 trajectories each, one with the TD-BA model, and other two with exact nonadiabatic coupling vectors using h-adjusted and p-adjusted velocities.
For both molecules, the TD-BA model was applied with analytical second time derivatives from quadratic regression (ΔT = 0.4 fs) and δη = 0.1 au. The conditions associated to δε and δϖ parameters were not applied. See Subsection VI.A for the definition and discussion of these four parameters.
The initial conditions were sampled from a harmonic oscillator Wigner distribution of the nuclei. For ethylene, they were restricted to the 9.30 ± 0.25 eV excitation window39. For fulvene, the initial conditions were restricted to the 4.00 ± 0.34 eV window.
The MCSCF calculations were done with Columbus (version 7, 09-Oct-2020)41,42. Dynamics was done with Newton-X (version 2.2 build 15)7,43 interfaced to Columbus. DC-FSSH based on TD-BA couplings with all options discussed in this paper is available in Newton-X version 2.2 build 15 or higher, to be used with any of the electronic structure methods and programs interfaced to Newton-X.
All datasets for ethylene and fulvene, with all Newton-X input and output files, are freely available (see Data availability)44–46.
Additional TD-BA DC-FSSH calculations were also done with linear-response time-dependent density functional theory (TDDFT) for fulvene and thiophene. Both sets of simulations were performed with the 6-31G(d,p) basis set using the CAM-B3LYP functional47 for fulvene and ωB97XD functional48 for thiophene. The classical equations of motion were integrated with Δt = 0.1 fs timestep for fulvene (maximum of 60 fs) and Δt = 0.5 fs for thiophene (maximum of 300 fs). The quantum equations were integrated with Δτ = Δt/20, using interpolated electronic quantities between classical steps. The initial conditions were sampled from a harmonic oscillator Wigner distribution of the nuclei. For fulvene, they are the same used for complete active space self-consistent field (CASSCF). For thiophene, the initial conditions were restricted to the 6.2 ± 0.5 eV excitation window. The TD-BA model was applied with analytical second time derivatives from quadratic regression (ΔT = 0.4 fs) and δη = 0.1 au. The conditions associated to δε and εϖ parameters were not applied. Dynamics were done with Newton-X (version 2.2 build 15) interfaced to Turbomole 7.549,50. GAMESS is a free and open source program that could be used in place of Turbomole to do the TDDFT calculations51.
The statistical analysis followed the protocol proposed in Ref. 39 to compare datasets. For ethylene, mean values and error bars (95% confidence interval) for seven observables (adiabatic-population and oscillator strength time constants and five structural channel yields) were computed from the trajectories. For fulvene, mean values and error bars were computed for six observables (first decay time constant, S1 population at 20 fs, S1 population at 60 fs, and three structural yields). These quantities were used to calculate the overlap score λx for each observable and the mean overlap score Λ(1,2) for the average over all λx39. λx measures the probability that the observable predictions in two datasets overlap (for some confidence interval). In turn, Λ(1,2) is near one if the predictions of both datasets significantly overlap and zero if they do not.
VI. Implementation and tests
A. General aspects
The second time derivative in Equation (12) can be computed numerically. In our implementation, in the first two classical timesteps of the trajectory and after any hopping, the nonadiabatic coupling is assumed to be null. At the third step, the second time derivative is computed by finite differences with the backward O(Δt) approximation
From the fourth timestep on, we use the backward O(Δt2) approximation
Alternatively, the second time derivative can be evaluated analytically after carrying out a quadratic regression over a sequence of Np = ΔT/Δt points accumulated over a period ΔT of a trajectory integrated with timestep Δt. After the regression, the predicted energy gap as a function of time is ΔEJL ≈ at2 + bt + c and the second time derivative is
Again, the nonadiabatic coupling is assumed to be null in the first two timesteps of the trajectory and after a hopping, while data are accumulated to do the regression.
To reduce the instabilities in the derivatives caused by small energy discontinuities during dynamics (they are widespread in MCSCF propagation52), we adopted three "cleaning" conditions. Thus, in addition to the choice of second time derivative method (and the associated ΔT value in the case of using quadratic regression), these conditions introduce three parameters — δη, δε, and δϖ — whose values are discussed later in this section.
The first of such conditions checks the magnitude of the coupling variation in one timestep Δt. If the rate of this variation is larger than a parameter δη, the current coupling is discarded in favor of the coupling in the previous timestep:
The second condition eliminates large couplings at energy gaps that are too large. TD-BA couplings are computed only for energy gaps smaller than a threshold δε and assumed to be zero otherwise:
The third condition avoids large couplings arising from points with large energy gaps but an anomalously large second time-derivative. The coupling is assumed to be null if the product is larger than an δϖ parameter:
Before doing any dynamics with the TD-BA couplings, we investigated the effect of using either method for evaluating second time derivatives and adopting the three cleaning conditions. To do so, we computed TD-BA couplings with different parameters setups for all points of the ethylene’s h-adjusted dataset39. The results are summarized in Figure 2. This figure reports the root mean square deviation (RMSD) between the TD-BA coupling σJL and the exact computed for the h-adjusted dataset
where N is the total number of timesteps in all trajectories.
Figure 2. Root mean square deviation (RMSD) between time-dependent Baeck-An (TD-BA) and exact nonadiabatic couplings computed for the ethylene’s h-adjusted dataset.
(a) RMSD is given for different values of δη with analytical (from quadratic regression over ΔT = 0.4 fs) and numerical (via finite differences) second time derivatives. The energy gap threshold (δε) is 12 eV and δϖ = 0.01 au. (b) RMSD is given for different quadratic regression periods (ΔT) for δη = 0.1 au, δε = 12 eV, and δϖ = 0.01 au. (c) RMSD is given for different energy gap threshold (δε) for δη = 0.1 au, ΔT = 0.4 fs, and δϖ = 0.01 au. (d) RMSD is given for different δϖ values for δη = 0.1 au, ΔT = 0.4 fs, and δε = 12 eV. Note the log scales in (a) and (d).
The RMSD in Figure 2 shows that couplings computed from analytical second derivatives from a quadratic regression are always more accurate than couplings computed from numerical second derivatives. The difference between the two procedures is particularly remarkable for δη larger than 1 au (Figure 2a). Nonetheless, the period ΔT over which the quadratic regression is done should be very short, smaller than 0.5 fs (Figure 2b). This means that for dynamics with large timesteps, using the quadratic regression may not be advantageous. The parameters δε (Figure 2c) and δϖ (Figure 2d) play only a minor role in the BA couplings' accuracy.
Density maps comparing the magnitudes of BA couplings to exact nonadiabatic couplings for all points of the h-adjusted dataset for ethylene are shown in Figure 3. As expected, most couplings are near zero. The largest couplings for each pair of states tend to fall on the diagonal, which implies a good agreement between TD-BA and the exact couplings. The order of magnitude is generally well represented, with the S0-S2 coupling magnitudes ten times smaller than those of the other two transitions. The TD-BA model tends to overestimate small-magnitude couplings, especially for the S0-S2 interactions. The model also shows some underestimation for large couplings in S0-S1 and S1-S2 interactions.
Figure 3. Density maps comparing the magnitudes of time-dependent Baeck-An (TD-BA) and exact nonadiabatic couplings for all points of each pair of states in the ethylene’s h-adjusted dataset.
TD-BA; smoothed over ΔT = 0.4 fs, and with δη = 1 au, δε = 12 eV, and δϖ = 0.01 au. A low-level contour is shown in each graph to guide the eye. The maps were computed with 672,812 data points.
Figure 4 shows results comparing TD-BA and exact couplings for a single trajectory of the h-adjusted dataset. In this case, after starting in S1, ethylene hops to S2 after 12.5 fs and stays there until it comes back to S1 at 48.1 fs. It finally returns to S0 at 93 fs (Figure 4a). The TD-BA model delivers a semiquantitative agreement with the exact results for the coupling magnitudes between all three pairs of states (Figure 4b, c, and d). The model correctly predicts the time where each peak occurs, as well as gives the right shape (intensity and width) of the coupling peaks. However, the TD-BA model shows spurious spikes that even the cleaning conditions discussed above could not get rid of.
Figure 4. Benchmark tests on one trajectory (TRAJ1) from the ethylene’s h-adjusted dataset.
(a) Potential energies as a function of time. Exact and time-dependent Baeck-An (TD-BA) coupling magnitudes during this trajectory are shown in (b) for S0-S1, (c) for S0-S2, and (d) for S1-S2.
B. Potential use of TD-BA in TDDFT dynamics
Linear-response TDDFT is the workhorse of excited-state calculations of large molecular systems53. Nevertheless, there are some critical obstacles to simulate nonadiabatic dynamics based on this method, especially when dealing with the internal conversion to the ground state. Because of its underlying approximations54, TDDFT not only fails in describing the state crossing between the ground and the first excited state, but also predicts the wrong dimensionality of the branching space around the crossing55. Such handicaps have not precluded cautious use of TDDFT for studying internal conversion to the ground state. Indeed, thanks to Send and Furche10, S1/S0 nonadiabatic coupling vectors for TDDFT are now a common feature of several quantum chemistry software programs. Our own research protocol has been to run TDDFT dynamics until the energy S1/S0 gap drops below some threshold (about 0.1 eV), stop the trajectory propagation there, and assume that point to be the hopping point34. Here, we were interested in testing whether the TD-BA model would allow going beyond this protocol, with actual hoppings to the ground state. To do so, we performed exploratory simulations with fulvene and thiophene.
The first challenge is to obtain an acceptable (TD) DFT description in regions close enough to the crossing region between the ground and excited state so that the hops between them can be computed. Therefore, we increased the size of the grid in the DFT calculations to the maximum value available in Turbomole 7.5 (grid size 7) as an attempt to avoid self-consistent field (SCF) convergence problems and negative excitation energies. Although this strategy successfully solved SCF convergence problems, our trajectories for fulvene were systematically terminated during the linear-response calculation of excitation energies due to singlet instabilities at small S0/S1 energy gaps. To further avoid singlet instabilities, we used the Tamm-Dancoff approximation (TDA)56, as proposed by Tapavicza et al.57 We computed 25 trajectories starting in the S1 state, of which 22 did not show any hopping between S1 and S0. In the remaining three trajectories, the hopping occurred at energy gaps of 0.41, 2.22, and 1.65 eV. This set of results suggests that the coupling between S1 and S0 is not well described for fulvene under the TDA approximation. Moreover, the S1/S0 energy gaps as a function of time were systematically too large (compared to those at CASSCF), indicating that the TDA description of the branching space was inadequate. Therefore, neither TDDFT nor TDA-DFT were suitable to be used for fulvene dynamics.
To test how system-dependent this failure was, we turned to thiophene. Surface hopping dynamics of thiophene with TDA was previously reported by Fazzi et al.,58 who showed that TDA was able to describe the nonradiative relaxation process bringing thiophene to the S1/S0 crossing region. We computed 12 trajectories starting from the S2 state. From those, six hopped to the ground state with reasonably small energy gaps. Nevertheless, five prematurely ended due to singlet instability at the crossing region. We tried to avoid the singlet instability by decreasing the time step to 0.1 fs, hoping to have a hopping before reaching the instability point. However, the instability persisted.
Therefore, although the TD-BA model delivers an alternative way to compute the nonadiabatic couplings between the ground and excited states, its application to nonadiabatic dynamics at TDDFT level is discouraging since it still depends on the good description of the relevant states close to the crossing region. We expect similar negative results for dynamics based on algebraic diagrammatic construction to second-order (ADC(2))59, which tends to yield negative excitations near the S1/S0 crossing. Note, however, that TD-BA can still be helpful for both TDDFT and ADC(2) dynamics involving only internal conversion between excited states.
VII. Dynamics results
A. Ethylene
Ethylene dynamics has been extensively investigated with many theoretical approaches, including wavepacket dynamics60, on-the-fly quasi-diabatic dynamics61, multiple spawning62–65 multiconfigurational Ehrenfest66, and surface hopping32,67–72. The ultrafast nonadiabatic dynamics of ethylene is driven by torsional and pyramidalization modes until the molecule finds the S1/S0 crossing seam.
The adiabatic-state population evolution based on TD-BA couplings is shown in
Figure 5 (top). The S1 state is quickly depopulated toward S0 within 100 fs. In the first 20 fs, S1 also receives some of the population, but it quickly returns to the lower states. The S1 lifetime fitted with an exponential decay is 107 ± 9 fs, where the margin of error is computed for a 95% confidence interval.
Population evolutions with TD-BA and exact couplings are also compared in
Figure 5 (top). The results are encouraging, delivering a semiquantitative agreement between the two datasets. The S1 lifetime in the h-adjusted dataset is 114 ± 10 fs, meaning that the results agree within the error bars. This level of variation is expected even within different sets of FSSH with exact couplings. For example, the lifetime of FSSH dynamics with exact couplings but adjusting velocities in the p direction is 120 ± 11 fs. Nevertheless, there are quantitative differences between TD-BA and the h-adjusted dataset. Although the TD-BA dynamics perfectly describes the S2 population evolution, the S1 and S0 evolutions bear some significant differences between the two datasets, with the S1 population after 50 fs always lower in the h-adjusted dataset.
Figure 5. Time evolution of time-dependent Baeck-An (TD-BA) decoherence-corrected fewest switches surface hopping (DC-FSSH) dynamics of ethylene.
(Top) mean adiabatic populations of S0, S1, and S2; (bottom) mean oscillator strength between the current and the ground states. "Exact" corresponds to DC-FSSH dynamics with exact couplings from the h-adjusted dataset.
The current state's mean oscillator strength is shown in
Figure 5 (bottom) for TD-BA and the h-adjusted dataset. The time constant obtained from an exponential fitting is 14 ± 2 fs, in excellent agreement with the reference value, 15 ± 2 fs. TD-BA correctly describes the damped oscillations in ethylene dynamics, which have been previously predicted by wavepacket60 and surface hopping67,68 simulations, and experimentally measured73,74.
One of the most significant differences between the TD-BA and the h-adjusted datasets is the fraction of molecules still in the excited states after 200 fs (Table 1). While the dynamics with exact couplings put this quantity at 9%, TD-BA dynamics has almost twice this amount, 16%. This disagreement is out of the margin of error of ±3%. Beyond this divergence, the ground-state-product' yields show much better agreement. At the end of the simulations, dynamics with TD-BA couplings predicts 38% of ground-state CH2CH2 structures, while dynamics with exact couplings, 39%. According to the TD-BA and reference sets, the expected amounts of ethylidene (CH3CH) are 14% and 18%, respectively. Single-H dissociation yields are 28% and 31%. Double-H dissociation yields are 1% and 3%. All these figures agree within their margin of errors.
Table 1. Mean value and error bars (95% confidence interval) of time constants, populations, and structural yields, computed for ethylene and fulvene dynamics.
Ethylene’s results for the h-adjusted (reference) and p-adjusted datasets of Ref. 39 are shown as well. and are the overlap scores for each observable computed between time-dependent Baeck-An (TD-BA) and the respective dataset.
| Observable | TD-BA | h dataset | () | p dataset | () |
|---|
| Ethylene | | | | | |
| Population time τp (fs) | 107±9 | 114±10 | 0.50 | 120±11 | 0.07 |
| Osc. strength time τo (fs) | 14±2 | 15±2 | 0.73 | 17±2 | 0.02 |
| S2, S1 CH2CH2 yield (%) | 16±3 | 9±3 | 0.00 | 12±3 | 0.05 |
| S0 CH2CH2 yield (%) | 39±4 | 38±4 | 0.84 | 41±4 | 0.83 |
| S0 CH3CH yield (%) | 14±3 | 18±3 | 0.13 | 15±3 | 0.94 |
| S0 C2H3 + H yield (%) | 28±4 | 31±4 | 0.60 | 27±4 | 0.83 |
| S0 C2H2 + 2H yield (%) | 1±1 | 3±2 | 0.22 | 3±1 | 0.01 |
| Fulvene | | | | | |
| 1st decay time τ1 (fs) | 11±1 | 9±1 | 0.00 | 9±1 | 0.00 |
| S1 population at 20 fs (%) | 46±7 | 34±6 | 0.00 | 42±7 | 0.66 |
| S1 population at 60 fs (%) | 12±5 | 3±3 | 0.00 | 6±3 | 0.02 |
| Planar S1→S0 hop (%) | 85±5 | 87±5 | 0.81 | 84±5 | 0.93 |
| Tw-Stretched S1→S0 hop (%) | 5±3 | 6±3 | 0.86 | 8±4 | 0.38 |
| Tw-Shrunk S1→S0 hop (%) | 9±4 | 7±4 | 0.73 | 8±4 | 0.91 |
The hopping counting and the mean energy gap at the hopping time for all transitions are shown in Table 2 for the three datasets. Concerning the mean gap, TD-BA is in excellent agreement with the other two sets for all transitions excepting S0→S2. This discrepancy is not statistically relevant given that there are only a few of such rare transitions in each set. The number of hoppings between S1 and S0 (both directions) also agrees well between TD-BA and the other two sets. Nevertheless, the number of hoppings between S1 and S2 (both directions) is twice as large in TD-BA as in the other sets. With exact couplings, 43% of trajectories populated S2. These trajectories returned to S1 afterward and, in general, remained there until they hopped to S0. With TD-BA couplings, the fraction of trajectories populating S2 was higher, 63%. These trajectories usually hopped forth and back between S1 ad S2 twice before converting to S0.
Table 2. Number of hoppings (count) and hopping energy-gap absolute mean value and standard deviation for the three sets of trajectories of ethylene and fulvene.
| TD-BA set | h set | p set |
|---|
| Count | Mean (eV) | St. Dev. (eV) | Count | Mean (eV) | St. Dev. (eV) | Count | Mean (eV) | St. Dev. (eV) |
|---|
| Ethylene | | | | | | | | | |
| S2→S1 | 520 | 0.35 | 0.22 | 259 | 0.39 | 0.28 | 252 | 0.43 | 0.32 |
| S2→S0 | 38 | 4.63 | 0.42 | 18 | 4.93 | 0.47 | 20 | 4.91 | 0.97 |
| S1→S2 | 563 | 0.61 | 0.45 | 281 | 0.63 | 0.43 | 285 | 0.62 | 0.42 |
| S1→S0 | 546 | 0.51 | 0.63 | 618 | 0.65 | 0.88 | 613 | 0.61 | 0.85 |
| S0→S2 | 4 | 3.91 | 0.54 | 2 | 2.12 | 0.29 | 6 | 1.88 | 0.19 |
| S0→S1 | 160 | 0.71 | 0.81 | 181 | 0.76 | 0.86 | 187 | 0.70 | 0.82 |
| Fulvene | | | | | | | | | |
| S1→S0 | 273 | 0.25 | 0.27 | 273 | 0.28 | 0.28 | 278 | 0.25 | 0.25 |
| S0→S1 | 92 | 0.57 | 0.56 | 82 | 0.56 | 0.40 | 91 | 0.47 | 0.32 |
We can understand the excess of S1-S2 hoppings in the TD-BA approach by analyzing a single trajectory computed with exact and TD-BA couplings, as illustrated in
Figure 6. The comparison is made to a p-adjusted dataset trajectory, limiting the differences to the coupling itself, without the influence of the velocity adjustment after hopping.
Figure 6. Characterization of the first 50 fs of a single trajectory of ethylene.
Potential energies of the S1 and S2 states as a function of time with (a) exact p-adjusted dataset and (b) time-dependent Baeck-An (TD-BA) couplings. The dots indicate the current state every 0.5 fs. Violet dots indicate ππ* diabatic character, and green dots indicate π*(2) character. The absolute value of the exact and TD-BA v.h projection computed for the reference trajectory is shown in (c). The potential energy profile of S1 and S2 along the CC distance, keeping the other coordinates at the S0 minimum, is given in (d).
In the initial stage of the S1 dynamics, the CC stretching strongly stabilizes the zwitterionic state S2 (π*(2)). As shown in
Figure 6a, the S1 and S2 states approach each other every approximately 15 fs. S2 destabilizes again when the CC bond shrinks due to the oscillatory motion. This initial oscillation of the S1-S2 gap is a common feature of all trajectories. When dynamics is propagated with TD-BA couplings, the hopping probability from S1 (ππ*) to S2 (π*(2)) at the stretched structure is larger than when propagating with exact couplings. Thus, the TD-BA trajectories hop to S2 about twice more often than the reference trajectories.
The reason for the larger probability is connected to the surface topography. With the exact coupling, the projection h on v is modulated by the velocity magnitude. When the CC stretching is maximum, the velocity in this internal coordinate is null, and v.h is null as well, causing the zero at 15 fs that we can see in the orange curve (p-adjusted dataset) in
Figure 6c. TD-BA coupling fails to capture this dependence on |v|, yielding a larger coupling peak at the same time. Nevertheless, this discrepancy in the probability does not play a significant role in the dynamics — at least not for ethylene. After populating S2, the molecule quickly returns to S1 within 5 fs or less, and the population evolution is not affected, as we have seen in Figure 5.
Table 1 gives the overlap score for all time constants and structural yields discussed above. When the TD-BA dataset is compared to the h-adjusted dataset, the scores () are excellent for the S0 CH2CH2 yield (0.8) and oscillator strength time constant (0.7). For the population time constant (0.5) and single-H dissociation (0.6), the scores are fine. They are poor for the less populated channels, ethylidene (0.1) and double-H dissociation (0.2). The results do not agree at all for excited-state CH2CH2 (0.0).
When the TD-BA dataset is compared to the p-adjusted dataset , the lifetime agreement deteriorates (0.07 score for population and 0.02 for oscillator strength), but most of the structural figures improve (Table 1). There is an almost perfect agreement for ground-state CH2CH2, CH3CH, and single H dissociation yields between the two sets, with overlap scores above 0.8. This means that part of the divergence to the h-adjusted dataset is caused not by the TD-BA couplings but by the p-adjustment of the velocities, which directly impacts the fragmentation process.
The statistical analysis in this paper follows the methods and protocol defined proposed in Ref. 39 Thus, the current results can be directly compared to those previously published datasets. Table 3 compiles a list of mean overlap scores Λ(1,2), comparing different FSSH setups to the h-adjusted reference dataset. The idea is that the implementation of new methods and testing of different features can always be analyzed in the same way to build a hierarchical comparison of procedures. We see that TD-BA yields Λ(1,2) of 0.43 at the same level as that of the p-adjusted dataset (0.41). This is equivalent to say that the results of ethylene DC-FSSH dynamics with TD-BA couplings and with exact couplings adjusting the velocity in the p direction are statistically equivalent for 500 trajectories and 95% confidence interval. Nevertheless, this statement does not mean that TD-BA results should be accepted without questioning. The average overlap score varies between 0 (complete disagreement with reference data) and 1 (perfect agreement). An average overlap score of 0.43 means that significant differences may be expected in comparison to the reference data. However, these differences are not worse than when we use exact couplings and do velocity adjustment in the p direction. As we discussed in the previous paper39, if we reduce the number of trajectories in the ensemble, the margins of error grow, and all overlap scores improve. For a small but typical ensemble of 100 trajectories, the mean overlap score for TD-BA dynamics is about 0.7, which is an acceptable agreement with the reference data. Therefore, TD-BA dynamics seems to be a promising approach for exploratory dynamics with few trajectories.
Table 3. Mean overlap score Λ(1,2) computed for different fewest switches surface hopping (FSSH) setups.
Ethylene: 500 trajectories in each dataset and 95% confidence interval averaged over seven observables. Fulvene: 200 trajectories in each dataset and 95% confidence interval averaged over six observables. The h- and p-adjusted datasets of ethylene are from Ref. 39.
| Dataset | Λ(1,2) |
|---|
| Ethylene | Fulvene |
|---|
| h-adjusted (reference) | 1 | 1 |
| p-adjusted | 0.41 | 0.59 |
| TD-BA | 0.43 | 0.40 |
B. Fulvene
Fulvene nonadiabatic dynamics has been discussed before based on different simulation approaches, including wavepacket propagation in reduced-dimensionality75–79, full-dimensionality direct-surface wavepacket propagation80,81, multiple spawning32, and surface hopping32,82. Fulvene has an extended S1/S0 crossing seam from planar to 90° twisted structures36,82,83. The central feature of its photophysics is the ultrafast sub-100 fs decay to the ground state, with periodic returns to S1 and internal conversion spread over the entire crossing seam.
For the analysis that follows, it is helpful to define the mean torsional angle as the average of the absolute values of the four dihedral angles around the C–CH2 bond, which in degrees is:
To compare the results of different datasets, we have defined three regions in the C–CH2 torsion and distance space. The first region (Planar) includes near planar structures with a mean torsional angle smaller than 30°. The second region (Twisted-stretched) includes twisted structures (≥ 30°) with a C–CH2 distance larger than 1.55 Å. The third region (Twisted-shrunk) includes twisted structures with C–CH2 distances smaller than 1.55 Å. The count of S1→S0 hoppings in each region is given in Table 1. These three regions are indicated in Figure 7a-c.
Figure 7. Characterization of fulvene dynamics.
C–CH2 mean dihedral angle and the CC distance at the S1→S0 hopping for the (a) time-dependent Baeck-An (TD-BA), (b) h-adjusted, and (c) p-adjusted datasets. The colors indicate the hopping time. The dotted circle shows the Franck-Condon (FC) position, and the crosses show the positions of the planar conical intersection (0°), the minimum on the crossing seam (63°), and the twisted conical intersection (90°). The dashed curve indicates the minimum on the crossing seam optimized with constrained angular values. The vertical and horizontal lines split the space into three regions, Planar (ϕC–CH2 < 30°), Twisted-stretched (ϕC–CH2 ≥ 30°, dC–CH2 ≥ 1.55 Å), and Twisted-shrunk (ϕC–CH2 ≥ 30°, dC–CH2 < 1.55 Å). Panel (d) shows the S1 population as a function of time for the three datasets.
The three datasets show that same general behavior. After the excitation, fulvene returns to the ground state within 10 fs, with S1 repopulations (recurrences) at 20, 40, and 60 fs (Figure 7d). Internal conversion to S0 can occur at different points of the crossing seam. As shown in
Figure 7 (a-c), the first wave of S1→S0 hoppings occurs in the Planar region, in the first 20 fs after excitation. In the h-adjusted dataset (Figure 7b), the mean value of this angle is 13° (±7°), and the mean C–CH2 distance is 1.58 Å (±0.06 Å). For comparison, the crossing seam has a maximum at 0° and 1.58 Å, implying that the first wave of S1→S0 hoppings occurs near this maximum. The second S1→S0 hopping wave happens between 20 and 40 fs at stretched C-CH2 distances (> 1.55 Å). In this second wave, the structure may be either planar or twisted. A third S1→S0 hopping wave happens between 40 and 60 fs, mostly at the Twisted-shrunk region. This third wave occurs near the minimum of the crossing seam, located at 63° and 1.48 Å.
To proceed with the statistical analysis, including the calculation of overlap scores, the S1 population error bar (95% confidence interval) was estimated with bootstrap through 10,000 repetitions of the 200 trajectories of each dataset.
The overlap scores are above 0.7 for the number of S1→S0 hoppings in each of the three regions of the torsion × distance space (Table 1). It means that the TD-BA model quantitatively agrees with both h- and p-adjusted datasets for predicting these observables. The only exception is the divergence in the number of Twisted-stretched hoppings between TD-BA and the p-adjusted dataset (0.38 overlap score). The overlap scores for the time and population observables do not agree so well, with values zero or near zero. Nevertheless, this divergence in the population evolution is not critical for the validity of the TD-BA model. As shown in Figure 8, the TD-BA dynamics follow the same qualitative trends in the h- and p-adjusted datasets. However, the initial TD-BA evolution is slightly slower, 11 ± 1 fs against 9 ± 1 fs of other datasets. Moreover, the TD-BA S1 population tends to be higher, especially when compared to the h-adjusted dataset. When considering all six observables, the mean overlap score Λ(1,2) between TD-BA and the h-adjusted reference dataset is 0.40 (Table 3), not much worse than the 0.59 obtained in the comparison between the p-adjusted dataset and the h-adjusted dataset.
Figure 8. Fulvene S1 population.
Top: time-dependent Baeck-An (TD-BA) and h-adjusted dataset; bottom: TD-BA and p-adjusted dataset. In each case, the shaded area indicates the 95% confidence interval sampled with bootstrap through 10,000 repetitions of the 200 trajectories of each dataset.
As a final remark, in Ref. 32, Ibele and Curchod reported surface hopping dynamics results for fulvene, also computed at CAS(6,6) level. Their S1 population at 20 fs computed with p-adjusted DC-FSSH was, like ours, near 0.4 (see
Figure 4 of that paper). Nevertheless, with h-adjusted DC-FSSH, this observable was reduced to about 0.2, which disagrees with our results. Although we cannot provide a definitive explanation for the divergence, it may be that their trajectory ensemble composed of 18 initial conditions repeated 10 times each is too biased to deliver a reliable population.
VIII. Conclusions
In Ref. 27, Baeck and An showed that the first-order nonadiabatic coupling vector could be approximated by a 1D model, depending only on the energy gap and the second derivative of this gap with respect to the coupling coordinate (see Equation (1)). In this paper, the BA approximation is recast on the time domain, giving rise to the TD-BA approximation (Equation (12)), which can be applied to multidimensional systems. The TD-BA couplings are equivalent to a coupling projection on the nuclear velocities, making them particularly suited to be employed in DC-FSSH simulations. TD-BA couplings immediately enable nonadiabatic dynamics with any electronic method that can provide excitation energies and energy gradients.
We have used the methodological protocols proposed in Ref. 39 to evaluate dynamics based on TD-BA couplings. Extensive simulations for ethylene and fulvene at MCSCF level showed that TD-BA DC-FSSH provides a correct qualitative picture of the dynamics and delivered results in semiquantitative agreement with the reference datasets computed with exact nonadiabatic couplings. The results' quality is statistically not much worse than performing DC-FSSH adjusting velocities after hopping in the momentum direction. Nevertheless, the model fails to adequately describe the nonadiabatic dynamics in regions with strong nonadiabatic couplings but small velocities (which are themselves challenging regions for surface hopping in general84).
In our implementation of the DC-FSSH with TD-BA couplings, we tested two ways of computing second time derivatives of the energy gaps and proposed three conditions to be satisfied to reduce numerical instabilities in the calculations. Although this initial modeling has delivered adequate results, there is a margin for further improvements, especially in detecting artifacts caused by small discontinuities in the potential energy surface.
We have also tested the TD-BA DC-FSSH with TDDFT to check whether this method would allow predicting hoppings to the ground state. In this case, however, the answer was negative due to singlet instabilities. This problem is not restricted to TD-BA and should happen in TDDFT dynamics based on exact couplings too. Dynamics based on TD-BA with TDDFT may still help describe internal conversion between excited states.
Given the uncertainties introduced by the TD-BA approximation, DC-FSSH dynamics with these couplings should be restricted to exploratory dynamics, where high accuracy is not a strong requirement. In this work, we have tested the model for the nonadiabatic dynamics of ethylene and fulvene. Although together these two molecules are representative of many typical excited-state cases, dynamics with TD-BA couplings should be thoroughly tested before being applied to other types of topography and densities, like in superexchange85, spin-exchange86, solvent-solute charge-transfer87, three-states crossing88, and dissociative89 internal conversion.
Faculty Opinions recommendedReferences
- 1.
Nelson TR, White AJ, Bjorgaard JA, et al.:
Non-adiabatic Excited-State Molecular Dynamics: Theory and Applications for Modeling Photophysics in Extended Molecular Materials.
Chem Rev.
2020; 120(4): 2215–2287. PubMed Abstract
| Publisher Full Text
- 2.
Crespo-Otero R, Barbatti M:
Recent Advances and Perspectives on Nonadiabatic Mixed Quantum-Classical Dynamics.
Chem Rev.
2018; 118(15): 7026–7068. PubMed Abstract
| Publisher Full Text
- 3.
Malone W, Nebgen B, White A, et al.:
NEXMD Software Package for Non-adiabatic Excited State Molecular Dynamics Simulations.
J Chem Theory Comput.
2020; 16(9): 5771–5783. PubMed Abstract
| Publisher Full Text
- 4.
Weingart O, Nenov A, Altoè P, et al.:
COBRAMM 2.0 - A software interface for tailoring molecular electronic structure calculations and running nanoscale (QM/MM) simulations.
J Mol Model.
2018; 24(9): 271. PubMed Abstract
| Publisher Full Text
- 5.
Humeniuk A, Mitric R:
DFTBaby: A Software Package for Non-adiabatic Molecular Dynamics Simulations Based on Long-Range Corrected Tight-binding TD-DFT(B).
Comput Phys Commun.
2017; 221: 174–202. Publisher Full Text
- 6.
Du L, Lan Z:
An On-the-Fly Surface-Hopping Program JADE for Nonadiabatic Molecular Dynamics of Polyatomic Systems: Implementation and Applications.
J Chem Theory Comput.
2015; 11(4): 1360–1374. PubMed Abstract
| Publisher Full Text
- 7.
Barbatti M, Ruckenbauer M, Plasser F, et al.:
Newton-X: a surface-hopping program for nonadiabatic molecular dynamics.
WIREs: Comp Mol Sci.
2014; 4(1): 26–33. Publisher Full Text
- 8.
Akimov AV, Prezhdo OV:
The PYXAID Program for Non-Adiabatic Molecular Dynamics in Condensed Matter Systems.
J Chem Theory Comput.
2013; 9(11): 4959–4972. PubMed Abstract
| Publisher Full Text
- 9.
Richter M, Marquetand P, González-Vázquez J, et al.:
SHARC: Ab Initio Molecular Dynamics with Surface Hopping in the Adiabatic Representation Including Arbitrary Couplings.
J Chem Theory Comput.
2011; 7(5): 1253–1258. PubMed Abstract
| Publisher Full Text
- 10.
Send R, Furche F:
First-order nonadiabatic couplings from time-dependent hybrid density functional response theory: Consistent formalism, implementation, and performance.
J Chem Phys.
2010; 132(4): 044107. PubMed Abstract
| Publisher Full Text
- 11.
Dallos M, Lischka H, Shepard R, et al.:
Analytic Evaluation of Nonadiabatic Coupling Terms at the MR-CI Level. II. Minima on the Crossing Seam: Formaldehyde and the Photodimerization of Ethylene.
J Chem Phys.
2004; 120(16): 7330–7339. PubMed Abstract
| Publisher Full Text
- 12.
Park JW, Shiozaki T:
Analytical Derivative Coupling for Multistate CASPT2 Theory.
J Chem Theory Comput.
2017; 13(6): 2561–2570. PubMed Abstract
| Publisher Full Text
- 13.
Gonon B, Perveaux A, Gatti F, et al.:
On the applicability of a wavefunction-free, energy-based procedure for generating first-order non-adiabatic couplings around conical intersections.
J Chem Phys.
2017; 147(11): 114114. PubMed Abstract
| Publisher Full Text
- 14.
Köppel H, Schubert B:
The concept of regularized diabatic states for a general conical intersection.
Mol Phys.
2006; 104: 1069–1079. Publisher Full Text
- 15.
Yu L, Xu C, Lei Y, et al.:
Trajectory-based nonadiabatic molecular dynamics without calculating nonadiabatic coupling in the avoided crossing case: trans↔cis photoisomerization in azobenzene.
Phys Chem Chem Phys.
2014; 16(47): 25883–25895. PubMed Abstract
| Publisher Full Text
- 16.
Zhu C, Nobusada K, Nakamura H:
New implementation of the trajectory surface hopping method with use of the Zhu–Nakamura theory.
J Chem Phys.
2001; 115(7): 3031–3044. Publisher Full Text
- 17.
Suchan J, Janos J, Slavicek P:
Pragmatic Approach to Photodynamics: Mixed Landau-Zener Surface Hopping with Intersystem Crossing.
J Chem Theory Comput.
2020; 16(9): 5809–5820. PubMed Abstract
| Publisher Full Text
- 18.
Xie W, Sapunar M, Došlic N, et al.:
Assessing the performance of trajectory surface hopping methods: Ultrafast internal conversion in pyrazine.
J Chem Phys.
2019; 150(15): 154119. PubMed Abstract
| Publisher Full Text
- 19.
Belyaev AK, Lasser C, Trigila G:
Landau–Zener Type Surface Hopping Algorithms.
J Chem Phys.
2014; 140(22): 224108. PubMed Abstract
| Publisher Full Text
- 20.
Tully JC:
Molecular-Dynamics with Electronic-Transitions.
J Chem Phys.
1990; 93: 1061–1071. Publisher Full Text
- 21.
Hammes‐Schiffer S, Tully JC:
Proton Transfer in Solution: Molecular Dynamics with Quantum Transitions.
J Chem Phys.
1994; 101(6): 4657–4667. Publisher Full Text
- 22.
Plasser F, Ruckenbauer M, Mai S, et al.:
Efficient and Flexible Computation of Many-Electron Wave Function Overlaps.
J Chem Theory Comput.
2016; 12(3): 1207–1219. PubMed Abstract
| Publisher Full Text
| Free Full Text
- 23.
Ryabinkin IG, Nagesh J, Izmaylov AF:
Fast Numerical Evaluation of Time-Derivative Nonadiabatic Couplings for Mixed Quantum-Classical Methods.
J Phys Chem Lett.
2015; 6(21): 4200–4203. PubMed Abstract
| Publisher Full Text
- 24.
Meek GA, Levine BG:
Evaluation of the Time-Derivative Coupling for Accurate Electronic State Transition Probabilities from Numerical Simulations.
J Phys Chem Lett.
2014; 5(13): 2351–2356. PubMed Abstract
| Publisher Full Text
- 25.
Granucci G, Persico M, Toniolo A:
Direct Semiclassical Simulation of Photochemical Processes with Semiempirical Wave Functions.
J Chem Phys.
2001; 114(24): 10608–10615. Publisher Full Text
- 26.
Westermayr J, Gastegger M, Marquetand P:
Combining SchNet and SHARC: The SchNarc Machine Learning Approach for Excited-State Dynamics.
J Phys Chem Lett.
2020; 11(10): 3828–3834. PubMed Abstract
| Publisher Full Text
| Free Full Text
- 27.
Baeck KK, An H:
Practical approximation of the non-adiabatic coupling terms for same-symmetry interstate crossings by using adiabatic potential energies only.
J Chem Phys.
2017; 146(6): 064107. PubMed Abstract
| Publisher Full Text
- 28.
An H, Baeck KK:
Practical and reliable approximation of nonadiabatic coupling terms between triplet electronic states using only adiabatic potential energies.
Chem Phys Lett.
2018; 696: 100–105. Publisher Full Text
- 29.
Westermayr J, Faber FA, Christensen AS, et al.:
Neural networks and kernel ridge regression for excited states dynamics of CH2NH2+: From single-state to multi-state representations and multi-property machine learning models.
Mach Learn Sci Technol.
2020; 1(1): 025009. Publisher Full Text
- 30.
Hu D, Xie Y, Li X, et al.:
Inclusion of Machine Learning Kernel Ridge Regression Potential Energy Surfaces in On-the-Fly Nonadiabatic Molecular Dynamics Simulation.
J Phys Chem Lett.
2018; 9(11): 2725–2732. PubMed Abstract
| Publisher Full Text
- 31.
Dral PO, Barbatti M, Thiel W:
Nonadiabatic Excited-State Dynamics with Machine Learning.
J Phys Chem Lett.
2018; 9(19): 5660–5663. PubMed Abstract
| Publisher Full Text
| Free Full Text
- 32.
Ibele LM, Curchod BFE:
A molecular perspective on Tully models for nonadiabatic dynamics.
Phys Chem Chem Phys.
2020; 22(27): 15183–15196. PubMed Abstract
| Publisher Full Text
- 33.
Gozem S, Melaccio F, Valentini A, et al.:
Shape of Multireference, Equation-of-Motion Coupled-Cluster, and Density Functional Theory Potential Energy Surfaces at a Conical Intersection.
J Chem Theory Comput.
2014; 10(8): 3074–3084. PubMed Abstract
| Publisher Full Text
- 34.
Abiola TT, Rodrigues NDN, Ho C, et al.:
New Generation UV-A Filters: Understanding Their Photodynamics on a Human Skin Mimic.
J Phys Chem Lett.
2021; 12(1): 337–344. PubMed Abstract
| Publisher Full Text
- 35.
Serrano-Andrés L, Merchán M, Borin AC:
Adenine and 2-aminopurine: Paradigms of modern theoretical photochemistry.
Proc Natl Acad Sci U S A.
2006; 103(23): 8691–8696. PubMed Abstract
| Publisher Full Text
| Free Full Text
- 36.
Blancafort L:
Photochemistry and Photophysics at Extended Seams of Conical Intersection.
Chemphyschem.
2014; 15(15): 3166–3181. PubMed Abstract
| Publisher Full Text
- 37.
Köppel H, Gronki J, Mahapatra S:
Construction scheme for regularized diabatic states.
J Chem Phys.
2001; 115(6): 2377–2388. Publisher Full Text
- 38.
Granucci G, Persico M:
Critical appraisal of the fewest switches algorithm for surface hopping.
J Chem Phys.
2007; 126(13): 134114. PubMed Abstract
| Publisher Full Text
- 39.
Barbatti M:
Velocity Adjustment in Surface Hopping: Ethylene as a Case Study of the Maximum Error Caused by Direction Choice.
J Chem Theory Comput.
2021. PubMed Abstract
| Publisher Full Text
- 40.
Ditchfield R, Hehre WJ, Pople JA:
Self-Consistent Molecular-Orbital Methods. IX. an Extended Gaussian-Type Basis for Molecular-Orbital Studies of Organic Molecules.
J Chem Phys.
1971; 54(2): 724–728. Publisher Full Text
- 41.
Lischka H, Shepard R, Müller T, et al.:
The generality of the GUGA MRCI approach in COLUMBUS for treating complex quantum chemistry.
J Chem Phys.
2020; 152(13): 134110. PubMed Abstract
| Publisher Full Text
- 42.
Lischka H, Shepard R, Shavitt I, et al.:
COLUMBUS, an ab initio electronic structure program, release 7.0. 2017. Reference Source
- 43.
Barbatti M, Granucci G, Ruckenbauer M, et al.:
NEWTON-X: A package for Newtonian Dynamics Close to the Crossing Seam (v. 2.2). 2018.
- 44.
do Casal MT, Toldo JM, Pinheiro M, et al.:
Fulvene DC-FSSH.
figshare.
Dataset. 2021. http://www.doi.org/10.6084/m9.figshare.14446998.v1
- 45.
Barbatti M:
Ethylene DC-FSSH: TD-BA.
figshare.
Dataset. 2021. http://www.doi.org/10.6084/m9.figshare.14010311.v1
- 46.
Barbatti M:
Ethylene DC-FSSH: Momentum Adjustment.
figshare.
Dataset. 2021. http://www.doi.org/10.6084/m9.figshare.13522856.v2
- 47.
Yanai T, Tew DP, Handy NC:
A New Hybrid Exchange-Correlation Functional Using the Coulomb-Attenuating Method (CAM-B3LYP).
Chem Phys Lett.
2004; 393(1–3): 51–57. Publisher Full Text
- 48.
Chai JD, Head-Gordon M:
Long-Range Corrected Hybrid Density Functionals With Damped Atom-Atom Dispersion Corrections.
Phys Chem Chem Phys.
2008; 10(44): 6615–6620. PubMed Abstract
| Publisher Full Text
- 49.
Balasubramani SG, Chen GP, Coriani S, et al.:
TURBOMOLE: Modular program suite for ab initio quantum-chemical and condensed-matter simulations.
J Chem Phys.
2020; 152(18): 184107. PubMed Abstract
| Publisher Full Text
| Free Full Text
- 50.
TURBOMOLE V7.5 2020, a development of University of Karlsruhe and Forschungszentrum Karlsruhe GmbH, 1989-2007, TURBOMOLE GmbH, since. 2007.
- 51.
Barca GMJ, Bertoni C, Carrington L, et al.:
Recent developments in the general atomic and molecular electronic structure system.
J Chem Phys.
2020; 152(15): 154102. PubMed Abstract
| Publisher Full Text
- 52.
Dawes R, Jasper AW, Tao C, et al.:
Theoretical and Experimental Spectroscopy of the S2 State of CHF and CDF: Dynamically Weighted Multireference Configuration Interaction Calculations for High-Lying Electronic States.
J Phys Chem Lett.
2010; 1: 641–646. Publisher Full Text
- 53.
Huix-Rotllant M, Ferré N, Barbatti M:
Time-Dependent Density Functional Theory: from the Fundamentals to Nonadiabatic Dynamics. In Quantum Chemistry and Dynamics of Excited States: Methods and Applications, González, L Lindh, R., Eds. John Wiley & Sons. 2020; 15–46. Reference Source
- 54.
Casida ME, Huix-Rotllant M:
Progress in Time-Dependent Density-Functional Theory.
Annu Rev Phys Chem.
2012; 63: 287–323. PubMed Abstract
| Publisher Full Text
- 55.
Levine BG, Ko C, Quenneville J, et al.:
Conical intersections and double excitations in time-dependent density functional theory.
Mol Phys.
2006; 104(5–7): 1039–1051. Publisher Full Text
- 56.
Hirata S, Head-Gordon M:
Time-Dependent Density Functional Theory within the Tamm–dancoff Approximation.
Chem Phys Lett.
1999; 314(3–4): 291–299. Publisher Full Text
- 57.
Tapavicza E, Tavernelli I, Rothlisberger U, et al.:
Mixed Time-Dependent Density-Functional Theory/Classical Trajectory Surface Hopping Study of Oxirane Photochemistry.
J Chem Phys.
2008; 129(12): 124108. PubMed Abstract
| Publisher Full Text
- 58.
Fazzi D, Barbatti M, Thiel W:
Modeling Ultrafast Exciton Deactivation in Oligothiophenes via Nonadiabatic Dynamics.
Phys Chem Chem Phys.
2015; 17(12): 7787–7799. PubMed Abstract
| Publisher Full Text
- 59.
Dreuw A, Wormit M:
The Algebraic Diagrammatic Construction Scheme for the Polarization Propagator for the Calculation of Excited States.
WIREs: Comp Mol Sci.
2015; 5(1): 82–95. Publisher Full Text
- 60.
Viel A, Krawczyk RP, Manthe U, et al.:
Photoinduced dynamics of ethene in the N, V, and Z valence states: a six-dimensional nonadiabatic quantum dynamics investigation.
J Chem Phys.
2004; 120(23): 11000–11010. PubMed Abstract
| Publisher Full Text
- 61.
Zhou W, Mandal A, Huo P:
Quasi-Diabatic Scheme for Nonadiabatic On-the-Fly Simulations.
J Phys Chem Lett.
2019; 10(22): 7062–7070. PubMed Abstract
| Publisher Full Text
- 62.
Ben-Nun M, Quenneville J, Martínez TJ:
Ab Initio Multiple Spawning: Photochemistry from First Principles Quantum Molecular Dynamics.
J Phys Chem A.
2000; 104(22): 5161–5175. Publisher Full Text
- 63.
Mori T, Glover WJ, Schuurman MS, et al.:
Role of Rydberg States in the Photochemical Dynamics of Ethylene.
J Phys Chem A.
2012; 116(11): 2808–2818. PubMed Abstract
| Publisher Full Text
- 64.
Tao H, Allison TK, Wright TW, et al.:
Ultrafast internal conversion in ethylene. I. The excited state lifetime.
J Chem Phys.
2011; 134(24): 244306. PubMed Abstract
| Publisher Full Text
- 65.
Yu JK, Bannwarth C, Hohenstein EG, et al.:
Ab Initio Nonadiabatic Molecular Dynamics with Hole-Hole Tamm-Dancoff Approximated Density Functional Theory.
J Chem Theory Comput.
2020; 16(9): 5499–5511. PubMed Abstract
| Publisher Full Text
- 66.
Saita K, Shalashilin DV:
On-the-fly ab initio molecular dynamics with multiconfigurational Ehrenfest method.
J Chem Phys.
2012; 137(22): 22A506. PubMed Abstract
| Publisher Full Text
- 67.
Barbatti M, Ruckenbauer M, Lischka H:
The photodynamics of ethylene: a surface-hopping study on structural aspects.
J Chem Phys.
2005; 122(17): 174307. PubMed Abstract
| Publisher Full Text
- 68.
Barbatti M, Granucci G, Persico M, et al.:
Semiempirical molecular dynamics investigation of the excited state lifetime of ethylene.
Chem Phys Lett.
2005; 401(1–3): 276–281. Publisher Full Text
- 69.
Mališ M, Luber S:
Trajectory Surface Hopping Nonadiabatic Molecular Dynamics with Kohn-Sham ΔSCF for Condensed-Phase Systems.
J Chem Theory Comput.
2020; 16(7): 4071–4086. PubMed Abstract
| Publisher Full Text
- 70.
Pradhan E, Sato K, Akimov AV:
Non-adiabatic molecular dynamics with ΔSCF excited states.
J Phys Condens Matter.
2018; 30(48): 484002. PubMed Abstract
| Publisher Full Text
- 71.
Hollas D, Šištík L, Hohenstein EG, et al.:
Nonadiabatic Ab Initio Molecular Dynamics with the Floating Occupation Molecular Orbital-Complete Active Space Configuration Interaction Method.
J Chem Theory Comput.
2018; 14(1): 339–350. PubMed Abstract
| Publisher Full Text
- 72.
Sellner B, Barbatti M, Müller T, et al.:
Ultrafast Non-Adiabatic Dynamics of Ethylene Including Rydberg States.
Mol Phys.
2013; 111(16–17): 2439–2450. Publisher Full Text
- 73.
Kosma K, Trushin SA, Fuss W, et al.:
Ultrafast dynamics and coherent oscillations in ethylene and ethylene-d4 excited at 162 nm.
J Phys Chem A.
2008; 112(33): 7514–7529. PubMed Abstract
| Publisher Full Text
- 74.
Kobayashi T, Horio T, Suzuki T:
Ultrafast Deactivation of the ππ*(V) State of Ethylene Studied Using Sub-20 fs Time-Resolved Photoelectron Imaging.
J Phys Chem A.
2015; 119(36): 9518–9523. PubMed Abstract
| Publisher Full Text
- 75.
Grohmann T, Deeb O, Leibscher M:
Quantum separation of para- and ortho-fulvene with coherent light: The influence of the conical intersection.
Chem Phys.
2007; 338(2–3): 252–258. Publisher Full Text
- 76.
Belz S, Grohmann T, Leibscher M:
Quantum dynamical simulations for nuclear spin selective laser control of ortho- and para-fulvene.
J Chem Phys.
2009; 131(3): 034305. PubMed Abstract
| Publisher Full Text
- 77.
Alfalah S, Belz S, Deeb O, et al.:
Photoinduced quantum dynamics of ortho- and para-fulvene: Hindered photoisomerization due to mode selective fast radiationless decay via a conical intersection.
J Chem Phys.
2009; 130(12): 124318. PubMed Abstract
| Publisher Full Text
- 78.
Blancafort L, Gatti F, Meyer HD:
Quantum dynamics study of fulvene double bond photoisomerization: The role of intramolecular vibrational energy redistribution and excitation energy.
J Chem Phys.
2011; 135(13): 134303. PubMed Abstract
| Publisher Full Text
- 79.
Ruiz-Barragan S, Blancafort L:
Photophysics of fulvene under the non-resonant stark effect. Shaping the conical intersection seam.
Faraday Discuss.
2013; 163: 497–512. PubMed Abstract
| Publisher Full Text
- 80.
Mendive-Tapia D, Lasorne B, Worth GA, et al.:
Controlling the mechanism of fulvene S(1)/S(0) decay: switching off the stepwise population transfer.
Phys Chem Chem Phys.
2010; 12(48): 15725–15733. PubMed Abstract
| Publisher Full Text
- 81.
Mendive-Tapia D, Lasorne B, Worth GA, et al.:
Towards converging non-adiabatic direct dynamics calculations using frozen-width variational Gaussian product basis functions.
J Chem Phys.
2012; 137(22): 22A548. PubMed Abstract
| Publisher Full Text
- 82.
Bearpark MJ, Bernardi F, Olivucci M, et al.:
Can Fulvene S1 Decay Be Controlled? A CASSCF Study with MMVB Dynamics.
J Am Chem Soc.
1996; 118(22): 5254–5260. Publisher Full Text
- 83.
Bearpark MJ, Blancafort L, Paterson MJ:
Mapping the intersection space of the ground and first excited states of fulvene.
Mol Phys.
2006; 104(5–7): 1033–1038. Publisher Full Text
- 84.
Subotnik JE, Ouyang W, Landry BR:
Can We Derive Tully's Surface-Hopping Algorithm from the Semiclassical Quantum Liouville Equation? Almost, but Only with Decoherence.
J Chem Phys.
2013; 139(21): 214107. PubMed Abstract
| Publisher Full Text
- 85.
Wang L, Trivedi D, Prezhdo OV:
Global Flux Surface Hopping Approach for Mixed Quantum-Classical Dynamics.
J Chem Theory Comput.
2014; 10(9): 3598–3605. PubMed Abstract
| Publisher Full Text
- 86.
Bai S, Barbatti M:
Mechanism of Spin-Exchange Internal Conversion: Practical Proxies for Diabatic and Nonadiabatic Couplings.
J Chem Theory Comput.
2019; 15(3): 1503–1513. PubMed Abstract
| Publisher Full Text
- 87.
Barbatti M:
Photorelaxation Induced by Water-Chromophore Electron Transfer.
J Am Chem Soc.
2014; 136(29): 10246–10249. PubMed Abstract
| Publisher Full Text
- 88.
Kistler KA, Matsika S:
Three-state conical intersections in cytosine and pyrimidinone bases.
J Chem Phys.
2008; 128(21): 215102. PubMed Abstract
| Publisher Full Text
- 89.
Sobolewski AL, Domcke W, Dedonder-Lardeux C, et al.:
Excited-State Hydrogen Detachment and Hydrogen Transfer Driven by Repulsive 1πσ* States: A New Paradigm for Nonradiative Decay in Aromatic Biomolecules.
Phys Chem Chem Phys.
2002; 4(7): 1093–1100. Publisher Full Text
