The symmetrical quasi-classical approach to electronically nonadiabatic dynamics applied to ultrafast exciton migration processes in semiconducting polymers

Exciton dynamics, or excitation energy transfer (EET), in conjugated semiconducting polymers plays an important role in determining the efficiency of organic electronic materials.¹ While EET in these systems has conventionally been treated using Förster rate theory,^1,2 recent experiments for representative conjugated polymers have provided evidence for the quantum coherent nature of EET on ultrafast time scales.^3–5 Coherent EET typically involves the migration of spatially delocalized excitons extending across 2-10 monomer units.³ These “spectroscopic units”^6–8 are delimited by geometric defects that break the conjugation, such as kinks and torsional deformations.^9–12 In line with this picture, several experimental studies suggest that torsional fluctuations interfere with EET on a subpicosecond time scale.¹³ Additionally, the formation and relaxation of local exciton ground states (LEGS),^14–16 which are subject to exciton self-trapping,^14,17 are taken to lead to a contraction of the photogenerated excitons. However, the microscopic nature of the delocalized excitons and the detailed mechanism of EET and the related relaxation events on ultrafast time scales have remained controversial.^7,18

Against this background, efficient simulation techniques adapted to coherent EET are needed, which are able to account for nonadiabatic transitions in high-dimensional systems exhibiting strong electron-phonon coupling. Classical investigations of the Ehrenfest type¹⁵ and surface hopping variety¹⁹ have been undertaken. They are, however, of limited accuracy as far as the treatment of electron-phonon correlations is concerned. Full quantum dynamical studies using efficient multiconfigurational methods, i.e., the multiconfiguration time-dependent Hartree (MCTDH) method^20,21 and its multi-layer (ML-MCTDH) variant,^22–24 are accurate but have generally been restricted with regard to the system size.^25–28 For systems composed of up to 20-30 sites, recent ML-MCTDH simulations support experimental observations that the relaxation of geometric defects such as torsional modes drives the exciton transfer across monomer units¹⁷ and reveals the nature of the excitons in polymer materials to be polaronic.^17,26,28

Ultimately, though, the size and/or complexity of systems which may be treated with full quantum dynamical simulations is inherently limited by the exponential scaling of their computational cost with system size. Classical methods, on the other hand, are not so limited and may potentially scale linearly to treat much larger (and presumably more realistic) polymer models, perhaps consisting of hundreds of monomer units (or more). The question, of course—in addition to whether such efficiency is realizable—is whether such classical approaches are sufficiently accurate: here, whether they can reproduce the essential features of exciton dynamics captured in full quantum simulations of modest-sized polymer models. If so, then classical approaches could potentially serve as powerful tools for understanding the dynamical behavior of semiconducting polymers at a much larger scale, and perhaps even the variation in excitation transfer capability of one semiconducting polymer species relative to another, and under what conditions.

In the last several years, a new methodology^29–38 has been developed for classically simulating electronically nonadiabatic molecular processes which is based on combining a symmetrical quasiclassical (SQC) windowing procedure with the well-known classical vibronic (nuclear + electronic) Hamiltonian of Meyer and Miller (MM). In a sense, this SQC/MM approach may be viewed as a generalization of standard classical molecular dynamics (MD) simulation³⁷ to include the treatment of electronically nonadiabatic processes—i.e., those molecular processes which involve multiple electronic states/Born-Oppenheimer surfaces—and this view stems from the fact that in MM vibronic dynamics, the electronic quantum states are mapped to ordinary classical harmonic oscillator degrees of freedom (DOF) whose vibrational excitations represent the occupations of the various electronic states. Thus, in conjunction with the standard classical treatment of the nuclear DOF in MD simulations, the MM model presents a simple, fully consistent, classical description of the full vibronic system, with electronic + nuclear DOF being properly coupled together in a unified Hamiltonian dynamics framework.

Compared to a quantum mechanical (QM) description, the advantage of the classical MM model is, as stated, that it may potentially be applied to the treatment of very large/complex molecular systems. Nonetheless, in order to determine electronic quantum state (and coherence) information, it is still necessary to apply some sort of quantization model to the electronic DOF in MM dynamics, and it is this aspect which is addressed by the symmetrical quasiclassical (SQC) windowing procedure. Indeed, the SQC windowing procedure is found central to the accurate and efficient application of the classical MM dynamics model.

The core SQC concept²⁹ is straightforward: it is simply that initial and final trajectory variables should be “quantized” in the same manner—i.e., symmetrically with respect to time, in the spirit of microscopic reversibility. Surprisingly, this simple principle applied³⁰ to the electronic DOF of the MM Hamiltonian,³⁹ in combination with a corresponding matching of the zero-point energy (ZPE) in the electronic MM oscillators, leads to a model for nonadiabatic molecular processes that—though stunningly simple—has been demonstrated to be accurate in a variety of molecular application contexts. These include gas phase reactive scattering³⁰ and photochemistry models, condensed-phase spin-boson-type models of dissipative electronic dynamics,³⁰ electron transfer reactions,³¹ Frenkel-type site-exciton models of light-harvesting complexes,³⁴ etc. In all these application contexts, the SQC/MM methodology is seen capable of quantitatively describing electronic quantum coherence, from regimes of extremely strong³⁷ to extremely weak³⁶ coupling between the electronic states, and including the de-coherence of such effects as modulated by coupling to the nuclear DOF (and without resorting to any ad hoc decoherence correction factors as frequently encountered in surface hopping approaches).

Recent work on the SQC/MM approach has also included several methodological/algorithmic developments, specifically: (i) an extension³⁵ of the SQC model to extract the full electronic density matrix at no additional cost relative to a standard SQC/MM calculation; (ii) a new triangle SQC windowing scheme³⁶ which, though developed to deal with the weak-coupling limit, is found superior to the original square histogram SQC model in all coupling regimes; and, most recently, (iii) a new “kinematic momentum” representation³⁸ for adiabatic MM vibronic dynamics which eliminates—analytically, without approximation—the need to calculate second-derivative couplings (as in the adiabatic multi-state/channel Schrödinger equation). The kinematic momentum representation is employed in this work and it is advantageous to do so anytime an adiabatic electronic representation is desirable, the most obvious case being when ab initio Born-Oppenheimer potential energy surfaces (PESs) are generated “on-the-fly” in realistic simulations.

This paper describes the application of the SQC/MM approach to ultrafast exciton dynamics in oligothiophene (OT), which is a model system related to poly-(3-hexylthiophene) (P3HT), a frequently employed donor material in organic photovoltaics.¹ Current theoretical investigations of exciton dynamics in such species often rely on the use of model Hamiltonians, and in recent MCTDH and ML-MCTDH studies,^26–28 a Frenkel site-exciton model Hamiltonian for OT type systems was derived from electronic structure calculations for suitable oligomer species,⁴⁰ using the mapping approach introduced in Ref. 41. The present work treats this same model Hamiltonian via the SQC/MM procedure, initially to once again assess the capabilities of the SQC/MM methodology but also with an aim toward treating much larger and more complex models of semiconducting polymers. For these purposes, a new protocol for extracting diabatic electronic populations from adiabatic SQC/MM calculations is developed, which is used to leverage the fact that a typical nonadiabatic process only involves a few low-lying adiabatic electronic states, the remaining states not requiring explicit treatment.

For benchmarking, SQC/MM calculations are compared with exact quantum calculations obtained via the multi-layer multi-configuration time-dependent Hartree (ML-MCTDH) method^22–24 and also with results calculated using the traditional Ehrenfest method, which can be viewed as a “baseline” classical treatment previously applied to similar model systems.¹⁵ The full quantum benchmark results (obtained via ML-MCTDH) are closely analogous to the simulations of Ref. 17, which provide evidence for the dynamics of self-trapping and torsionally induced exciton migration in a low-temperature regime. Below, it is shown that the SQC/MM simulation results agree well with the benchmark ML-MCTDH calculations and significantly outperform the Ehrenfest simulations.

The remainder of the paper is organized as follows: In Sec. II, the model Frenkel/site-exciton OT Hamiltonian is presented, as well as the details of the ML-MCTDH, Ehrenfest, and SQC/MM methods used to treat it. The results of these simulations are shown and discussed in Sec. III. Section IV concludes this paper.

The exciton dynamics study presented in this paper centers on the treatment—by both quantum and classical mechanics—of a Frenkel/site-exciton Hamiltonian model corresponding to a 20 monomer/site OT segment which was previously parametrized via ab initio electronic structure calculations.^17,41 This OT Hamiltonian is described first, followed by a detailed description of the classical SQC/MM dynamics protocol used to treat it. The benchmark quantum calculations were performed by the ML-MCTDH^22–24 methodology, and relevant details of these calculations will be discussed as they arise.

Following Refs. 17, 26, and 28, the diabatic exciton model Hamiltonian for the OT segment is formulated on the basis of Frenkel configurations, which are taken to correspond to single-monomer excitations in the direct-product space spanned by F monomer units.

The total Hamiltonian is

i.e., the nuclear kinetic energy, given for generality in terms of the metric tensor G, plus the electronic Hamiltonian, given (in this case) by

which combines state/site-dependent potentials ${ V k site }$⁠, each corresponding to an electronic configuration localizing an excitation on the kth of the F monomers, a state-independent bath component V_bath, and excitonic (diabatic) couplings ${ V k k + 1 exc }$ here taken to link electronic states corresponding to excitations on adjacent monomers. Several types of nuclear modes are included, notably, a set of site-local ring-breathing modes, X ≡ {X_k}, along with two types of “correlated” modes that connect neighboring monomer sites: inter-ring bond-stretch modes, Y ≡ {Y_kk+1}, and torsional modes, θ ≡ {θ_kk+1}. These are schematically illustrated in Fig. 1. Note the dual indices of the “correlated” modes indicating to which two monomers each mode corresponds: i.e., each monomer potential is a function of one ring-breathing mode (X_k), two bond-stretches (Y_kk−1, Y_kk+1), and torsions (θ_kk−1, θ_kk+1) shared with adjacent monomers.

The state/site-dependent potentials in Eq. (2) are expressed in terms of the given nuclear modes (X, Y, θ) as

with the ground-state potential

and the difference potential

where c_E is a constant monomer excitation energy and $v G$ and $v E$ are the ground (G) and excited (E) state monomer potentials, respectively—i.e., the kth electronic state consists of the excited state PES of the kth monomer plus the ground state PES of the other monomers. As detailed in the supplementary material, these monomer potentials are the same as used in Ref. 17—assumed separable in the monomer coordinates and generated by point-wise mapping of ab initio generated thiophene oligomer PESs to the desired monomer potentials according to the analytical procedure described in Ref. 41. As an additional simplification, it is assumed here that the monomer potentials are harmonic in the high-frequency ring-breathing and bond-stretch modes (the torsional angular dependence being represented in a cosine series).

The excitonic coupling in Eq. (2) is assumed to be independent of the nuclear DOF, i.e.,

which is a reasonable approximation in view of the results reported in Ref. 41.

Turning to the V_bath term in Eq. (2), the bath modes B ≡ {B_k} are also taken to be harmonic and, generally, would be coupled to one or more of the torsional modes. In this work, however, analogous to the calculations reported in Ref. 17 and to keep the ML-MCTDH calculations to a reasonable computational effort, only the inter-monomer torsion at the middle of the OT segment, θ_{10 11}, is treated dynamically (which, as explained below, serves to break the segment into two 10-mer “spectroscopic units”). Accordingly, only θ_{10 11} is coupled to the harmonic bath, which is thus given by

where, as usual (and assuming mass-weighted coordinates), the bath couplings {c_j} are adapted to an Ohmic spectral density, $c j = ω j ( 2 / π ) η Δ ω$ with η being the friction coefficient and Δω being the frequency spacing that determines the Poincaré recurrence time T_P = 2π/Δω. In this work, 10 bath modes were employed (N_bath = 10). Additional details are given in the supplementary material.

Finally, it is noted that the explicit form of the nuclear kinetic energy $1 2 P ^ T G P ^$ is adopted from Refs. 17 and 42, with the elements of the metric tensor G given in terms of the masses and moments of inertia for the various nuclear modes (X, Y, θ, B). Its details are also provided in the supplementary material.

With the model OT system as just described, the goal is to simulate vibronic dynamics after the photoexcitation of a particular contiguous/conjugated subset of monomers—i.e., to simulate the exciton migration, relaxation, and self-trapping after the photoexcitation of a particular “spectroscopic unit” as defined by a defect in the OT polymer chain. In these simulations, as stated above, the spectroscopic unit-defining defect is chosen to be in the torsional angle between monomers 10 and 11. With an initial displacement of θ_{10 11} = 110°, the OT chain is put in a non-planar configuration, breaking the π-conjugation at θ_{10 11}, effectively defining two spectroscopic units as left and right 10-mer segments. Also, as stated above, θ_{10 11} is the only torsional bond angle treated dynamically; all others are locked at their ground-state equilibrium geometries of θ_eq = 160°.

Within the single-excitation Frenkel manifold, the lowest-lying electronic configurations corresponding to this initial nuclear geometry are nodeless excitons that are localized on one or the other side of the torsional defect. Here, this initial local exciton (LE) state is taken to be on the left segment (left-LE) and this, of course, must be properly set up for the dynamical simulations.

For the full quantum ML-MCTDH calculations, the initial left-LE electronic configuration is established via imaginary-time propagation in the electronic subspace, prior to the real-time quantum propagation.^17,26 This corresponds to the lowest adiabatic electronic eigenstate of the left 10-mer OT segment.

The electronic amplitudes emerging from the imaginary-time propagation are also used to set the initial conditions of the electronic variables in the Ehrenfest calculations. For the SQC/MM simulations, however, a more complicated—but nevertheless straightforward—procedure is required, as described in detail below.

The nuclear DOF are initialized in the electronic ground state, according to the ground-state potential V₀ of Eq. (4), except for the active torsional mode which is given an initial displacement at θ_{10 11} = 110°. The remaining bond angle torsions (θ) are fixed at their ground-state equilibrium geometry of θ_eq = 160°. Thus, the initial nuclear conditions correspond to a vertical (Frank-Condon) excitation into an electronic state with delocalization broken by the torsional defect.

The SQC/MM approach constitutes a very general methodology for classically simulating electronically nonadiabatic processes when what is sought is a quantitative treatment of the coupled vibronic (nuclear + electronic) dynamics, including electronic coherence/decoherence phenomena. The approach is thoroughly documented in a series of recent papers^29–38 with Ref. 37, in particular, providing a concise review of the methodology. A brief summary is now provided in terms of its two key components: (i) the classical Meyer-Miller (MM) Hamiltonian and (ii) the symmetrical quasi-classical (SQC) windowing model. The focus will be on the aspects of the SQC/MM approach relevant to the application at hand and, in particular, the application of appropriate initial and final electronic boundary conditions as dictated by the OT model and specified initial conditions.

The Meyer-Miller (MM) Hamiltonian³⁹ represents the dynamics of the electronic states in nonadiabatic processes as the classical motion of a set of “electronic oscillators” whose excitations represent the occupations of the various electronic states. In terms of diabatic potential energy surfaces (PESs) {H_kk(R)} associated with F electronic states—in general, dependent on nuclear coordinates R—and their off-diagonal couplings {H_kk′(R)}, it is given by

where {p_k, x_k} are the momenta and coordinates of the aforementioned F electronic oscillator DOF [their harmonic form is readily apparent in Eq. (8)] and γ is a zero-point energy (ZPE) parameter, originally taken to be $1 2$ (the quantum value) but, in the SQC approach, always taken to be $≈ 1 3$ which is found consistently optimal (see, e.g., Ref. 37). Note that the nuclear kinetic energy’s form in Eq. (8) is given in terms of the metric tensor G, as it is in Eq. (1), and again, see the supplementary material for its exact form.

Meyer and Miller originally proposed Eq. (8) as a classical Hamiltonian, a “classical analog” model; however, Stock and Thoss later showed⁴³ that if this Hamiltonian was treated quantum mechanically—i.e., as an operator in the time-dependent Schrödinger equation—it would actually give the exact quantum dynamics. In this sense, the classical dynamics generated from Eq. (8) constitutes more than just another model but in fact the well-defined classical-limit of the true multi-state vibronic quantum dynamics.

The time-dependent Schrödinger equation for multiple electronic states may also of course be written in an adiabatic representation—i.e., in terms of the Born-Oppenheimer PESs {E_k(R)} and first- and second-derivative coupling vectors between the different electronic states,

with {ϕ_k} being the adiabatic electronic wavefunctions/states—and, likewise, the MM Hamiltonian may also be written in an adiabatic representation,

where the derivative couplings {d_kk′} appear in the kinetic energy as a modification to the nuclear momentum as

Note that in Eq. (9), just as in the adiabatic Schrödinger equation, it is the derivative coupling vectors which are responsible for the electronic transitions (instead of off-diagonal elements of the potential matrix {H_kk′}, as in a diabatic representation).

Because the two forms of the MM Hamiltonian, diabatic [Eq. (8)] and adiabatic [Eq. (9)], are totally equivalent, one has the option of choosing whichever electronic representation/basis is most suitable (i.e., convenient and/or efficient) for the problem at hand. The Frenkel exciton Hamiltonian of Eq. (1) is given in a diabatic electronic representation having 20 electronic states, each corresponding to an excitation localized on one of the 20 monomers. However, the vibronic dynamics in these calculations, as is often the case (and as verified below), actually only involves a few of the lowest-lying adiabatic electronic states. Thus, the decision was made to take advantage of this common scenario and to perform simulations directly in the adiabatic representation. It is worth noting that the Frenkel exciton Hamiltonian for OT used here was parametrized via ab initio electronic structure calculations, which are necessarily performed in the adiabatic basis. It is also worth noting that by employing the adiabatic representation in these calculations, we leverage a very compact electronic basis and by working out how to do this effectively and consistently within the SQC/MM protocol, the ground work is laid for future simulations that could involve many more monomers and nuclear DOF.

One final important point regarding the use of the adiabatic representation concerns the possible presence of second-derivative coupling vectors {D_kk′}. These do not appear explicitly in the classical adiabatic Hamiltonian of Eq. (9), but they will nevertheless appear in the standard adiabatic equations of motion (EOM) from applying Hamilton’s equations to Eq. (9) (which is unsurprising as they also appear in the standard form of the adiabatic Schrödinger equation). This would likely be the most significant drawback of using the adiabatic MM Hamiltonian, for example, in an on-the-fly ab initio dynamics calculation, as second-derivative couplings are generally not available from standard ab initio quantum chemistry packages. Fortunately, however, as developed in Ref. 38, this drawback is totally eliminated—analytically, without approximation—by simply expressing the nuclear dynamics in a “kinematic momentum” representation. The details of this are given in  Appendix A including the full set of kinematic adiabatic EOM used exclusively in this work.

The classical MM Hamiltonians of Eqs. (8) and (9) (diabatic and adiabatic) do not by themselves imply any sort of quantization, and so this aspect of the procedure is provided by the symmetrical quasi-classical (SQC) windowing model, as now described.

The SQC quantization protocol operates in terms of classical “action” variables defined as

which are the classical/continuous analogs of the vibrational quantum numbers of the electronic oscillators in Eqs. (8) and (9). It is clear from Eqs. (8) and (9) that they represent electronic occupations—e.g., they are the pre-factors multiplying the diagonal PES matrix elements [H_kk(R) and E_k(R) in Eqs. (8) and (9), respectively]. Likewise, the classical “angle” variables canonically conjugate to the actions {n_k} are given by

and so the inverse transformation between the Cartesian oscillator variables {p_k, x_k} [in terms of which the actual dynamics are computed via Eqs. (8) and (9)] and action-angle oscillator variables {n_k, q_k} (in terms of which the windowing/quantization is done) is

The relationships in Eqs. (11) and (12) may also be expressed in the polar representation,

where the electronic amplitudes {c_k} have been identified. (See, e.g., Ref. 38.)

Typical applications of the SQC/MM approach have used the protocol to compute the time-dependent population dynamics of the various electronic states involved in the modeled electronically nonadiabatic molecular process—i.e., the relevant diagonal elements of the electronic density matrix {ρ_kk(t)}. Using the standard square histogram SQC window functions, an electronic quantum state k of, say, F electronic states is represented in terms of the set of F electronic actions n [as defined in Eq. (11a)] by

where $w$_N(n) windows a single electronic DOF/action n about the quantum value N and is defined as

That is, the SQC “quantization” model of Eqs. (14) and (15) defines quantum state k as requiring action variable n_k to be within ±γ of the quantum value n_k = 1 and all the other action variables {n_k′≠k} to be within ±γ of the quantum value n_k′ = 0. Furthermore, the parameter γ used in Eqs. (14) and (15) to define the width of the square histogram SQC windows is crucially taken to be the same parameter γ that defines the fractional ZPE in the MM Hamiltonians of Eqs. (8) and (9) [and in the transformations between action-angle and Cartesian variables in Eqs. (11) and (12)]; the value of γ = 0.366 is used exclusively in the present OT application, as has been done in all previous applications of the SQC/MM approach employing square histogram windowing functions (at least by these authors). (See the appendix of Ref. 30 for one theoretical justification of the significance of this particular value.)

As stated, Eq. (14) defines the diagonal elements of the electronic density matrix in the SQC model. More generally, it is also straightforward (and necessary in this work) to compute the off-diagonal electronic density matrix elements ρ_k≠k′(t). Following Ref. 35, this is done by representing these elements with an SQC histogram window function given by

with q as defined in Eq. (11b). Equation (16) thus manifests the explicit dependence of the electronic coherence between two states k and k′ on the difference in phase-angles q_k − q_k′ of the electronic oscillators in Eqs. (8) and (9), and note that Eq. (16) collects these phase differences when their conjugate actions fall within a region of the action space centered between two quantum states k, k′ at $( n k , n k ′ ) = ( 1 2 , 1 2 )$⁠.

Employing the window functions defined by Eqs. (14)–(16) and for convenience defining

calculation of F × F electronic density matrix elements {ρ_kk′(t)} at the time t propagated from an initial given electronic state i and initial nuclear distribution ρ (P₀, R₀) for the G nuclear DOF amounts to Monte Carlo evaluation of the phase-space average

for all k, k′ ∈ [1, F], followed by renormalization via

In other words, in an SQC/MM simulation, trajectories are

1. initialized via importance sampling: ρ (P₀, R₀) is used to sample the G nuclear DOF, and the F electronic actions are sampled according to W_i(n₀) with initial angles q₀ chosen randomly in the interval [0, 2π] and (n₀, q₀) converted to Cartesian (p₀, x₀) according to Eq. (12),

2. propagated to time t via Hamilton’s equations from the MM Hamiltonian of Eq. (8) (diabatically) or Eq. (9) (adiabatically), and finally

3. collected at time t via W_kk′(n, q) [Eq. (17)]—after conversion to action-angle variables (n_t, q_t) according to Eq. (11)—to calculate ${ ρ ̃ k k ′ ← i ( t ) }$ which are renormalized [Eq. (18b)] to {ρ_kk′←i(t)}.

Note that just one ensemble of trajectories may be used to calculate all matrix elements ${ ρ ̃ k k ′ ← i ( t ) }$ in Eq. (18) and also that the renormalization constant in Eq. (18b) only involves a sum over diagonal elements. This is the same procedure used in Ref. 35 and in all the other SQC/MM papers, Refs. 30–34 and 36–38, with the additional use of Eq. (16) to calculate the off-diagonal density matrix elements.

For the reasons described in Sec. II B 1, the decision was made in this work to exploit the very compact electronic representation made possible by the use of the adiabatic version of the MM Hamiltonian [Eq. (9)]. As a practical matter, there are only 2-4 adiabatic electronic states of significance (i.e., having any appreciable population) at any point during the simulation; therefore, one can perform adiabatic SQC/MM calculations here using just 1/5th the number of electronic states versus a site-local diabatic approach, with no appreciable sacrifice in accuracy.

However, despite the complete equivalence of the classical dynamics generated from either adiabatic or diabatic MM Hamiltonians—necessarily so, since one is a canonical transformation of the other—adiabatic and diabatic SQC boundary conditions are not the same—one must be converted to the other for the simulations to be equivalent—and so a unique aspect of the present work is the handling of different electronic basis/representations within the SQC windowing/quantization model. And, of course, this involves both initial and final boundary conditions: initially, electronic variables need to be set up in the left-LE electronic configuration, as described in Sec. II A 2, and at all final times t, diabatic electronic populations need to be extracted from the adiabatic dynamical calculations, this final step being because the goal here is to analyze exciton migration and localization, and the (Frenkel-exciton) diabatic representation describes this site-specific localization of electronic excitation.

As described in Sec. II A, the 20-mer OT system is to be initialized with a torsional defect at its center, splitting it into a pair of 10-mer “spectroscopic units,” whereby the initial electronic excitation is the (partially) localized exciton ground state of the left segment—i.e., the initial left-LE. The procedure for establishing this initial condition is as follows.

First, the nuclear coordinates and momenta are initialized from an appropriate distribution corresponding to the monomer ground state potentials. For the ring-breathing and bond-stretch modes (X, Y), two choices are made in this work: (i) nuclear coordinates and momenta are sampled from the Wigner distribution for the ground state (i.e., zero temperature, as chosen in the ML-MCTDH calculations) or (ii) the nuclear coordinates are initialized at exactly their minimum energy configuration with all nuclear momenta taken to be zero (i.e., δ-function distributions). Section III presents the results for both choices.

With the nuclear DOF initialized, the corresponding initial left-LE configuration can be established for the electronic DOF. Because the dynamics is to be run in the adiabatic representation, the protocol requires the nuclear coordinate-dependent adiabatic-to-diabatic transformation matrix obtained from diagonalizing the electronic Hamiltonian Ĥ_elec(R) [Eq. (2)]. This unitary transformation is designated

where ${ ϕ k adia ( R ) }$ are the adiabatic electronic basis functions parametrically dependent on the nuclear coordinates R ≡ {X, Y, θ, B} and, for simplicity, a bare index in a bra-ket refers to the underlying diabatic/Frenkel basis.

The left-LE initial condition, however, is not the ground adiabatic electronic state of the full Ĥ_elec but rather the lowest-lying electronic state with excitation localized on the left segment—i.e., it is the ground adiabatic state of the left 10-mer, which roughly corresponds to a coherent superposition of the two lowest-lying adiabatic states of the full 20-mer. This electronic state may be determined by block-diagonalizing Ĥ_elec in the diabatic basis with the coupling(s) between left- and right-segments set to zero,

The corresponding coordinate-dependent unitary matrix—transforming from this coordinate-dependent LE electronic basis ${ ϕ k L E }$ [the eigenvectors of Eq. (20)] to the diabatic/Frenkel basis—is designated as

The portion of the lowest-lying eigenstate of Eq. (20) which relates to excitation on the left segment is taken to be the initial electronic state sampled via the SQC window function of Eq. (14).

What remains is to convert the initially sampled left-LE electronic variables to the adiabatic representation for dynamics propagation. Having sampled the initial nuclear configuration from ρ (P₀, R₀) and determined the two transformation matrices of Eqs. (19) and (21), the initial sampled electronic variables (n₀, q₀) are converted to Cartesian (p₀, x₀) via Eq. (12) and then further transformed to the adiabatic basis by sequentially applying the foregoing unitary transforms,

the first transforming from left-LE to diabatic/Frenkel basis and the second (note the adjoint †) transforming from diabatic/Frenkel to adiabatic basis. Note that these linear operations are appropriate because the electronic amplitudes {c_k} are linearly related to the Cartesian oscillator variables by $c k = 1 2 ( x k + i p k )$ [Eq. (13)].

This completes the selection of the initial trajectory variables [step 1 in evaluating Eq. (18)] consistent with the left-LE boundary condition, after which propagation in the adiabatic representation may be performed [step 2 in evaluating Eq. (18)]. Because it only takes the coherent superposition of a few adiabatic eigenvectors to accurately represent the initial left-LE configuration and subsequent exciton dynamics—i.e., only a few of the adiabatic {c_k} are ever appreciably different from zero—only a few of the adiabatic electronic DOF need to be propagated in the MM dynamics and eventually windowed at each time t > 0 which greatly enhances the statistical convergence and efficiency of the approach.

The final step of the SQC protocol—after time-evolving the electronic and nuclear DOF—is to window the final electronic action-angle variables (n, q) [step 3 in evaluating Eq. (18)]. In the present work, a novel aspect of this last step stems from the fact that the chief interest is in exciton migration and localization over the 20-mer OT chain and, of course, this is naturally tied to the diabatic/Frenkel site-local electronic basis. Accordingly, some protocol must be devised for extracting site-local diabatic electronic populations from adiabatic MM dynamics.

Conceptually, the adiabatic density matrix ρ_kk′←i(t) is calculated at time t by evaluating Eq. (18) in the adiabatic basis with the initial density matrix taken to be a particular coherent superposition of a few adiabatic states which represent the left-LE. Furthermore, the electronic density matrix/operator, once computed, may in principle be changed from one basis to another by unitary transform (just like all QM matrices/operators). The difficulty, however, is that (in this case) the transformation matrix—the adiabatic-to-diabatic unitary matrix of Eq. (19)—is coordinate-dependent, meaning that a different transformation matrix applies to each trajectory at each time t.

Nevertheless, this scenario is handled consistently as follows: for each adiabatic trajectory, the window function W(n, q)—Eq. (17) in matrix form—is applied to the adiabatic actions and angles (n_t, q_t). This result, W(n_t, q_t), is what a particular trajectory would contribute to an adiabatic density matrix. The corresponding contribution to the diabatic density matrix is simply the same result, W(n_t, q_t), transformed to the diabatic basis using the adiabatic-to-diabatic transformation matrix of Eq. (19). That is, each trajectory propagated in the adiabatic basis makes a contribution to the diabatic/site-local density matrix given by

Summing the contributions of Eq. (23) from all trajectories and renormalizing yields the final desired diabatic density matrix. These steps are, of course, equivalent to replacing Eq. (18a) with

where the time-evolved actions and angles (n_t, q_t) on the RHS are adiabatic, but the LHS is the desired diabatic/Frenkel-exciton density matrix. Equations (23) and (18a′) show the need to calculate the full SQC electronic density matrix in these kinds of simulations.

One may view the strategy of Eq. (23), applied in Eq. (18a′), as a new general procedure within the SQC protocol for extracting diabatic electronic populations from adiabatic MM dynamics simulations. Interestingly, the technique is similar to what was done recently to extract diabatic populations from adiabatic surface hopping simulations,⁴⁴ and these methods also rely on computing off-diagonal elements of the electronic density matrix. In the context of the SQC/MM approach, however, the strategy of Eq. (23) is nevertheless untested and, therefore, to preliminarily validate it,  Appendix B presents its application to several standard spin-boson problems, where it is indeed shown that trajectory propagation and windowing in the adiabatic electronic basis, followed by extraction of diabatic populations via a coordinate-dependent adiabatic-to-diabatic transformation [Eq. (23)], does in fact produce results in very good agreement with those obtained via a standard SQC/MM diabatic calculation.

To recap, the foregoing diabatic extraction procedure in combination with the left-LE initialization described in Subsection II B 3 a may be summarized as

where, again, this expression replaces Eq. (18a) and is evaluated by Monte Carlo and renormalized in the usual manner described with respect to Eq. (18), except that the initial electronic DOF are sampled in the left-LE basis and transformed to the adiabatic basis with Eq. (22). Note that the upper limit of the summation, F′ = 4 (much less than F = 20), explicitly leverages the compact dimensionality of the relevant adiabatic electronic subspace.

Benchmark quantum dynamical calculations were performed using the ML-MCTDH method^23,24 (Heidelberg MCTDH package⁴⁵) in the diabatic representation of Eqs. (1) and (2) and comprising a seven-layer representation of the wavefunction, as detailed in the supplementary material. These calculations yield very similar results to those previously carried out in Ref. 17, the difference being that, in Ref. 17, more general anharmonic potentials were used for the ring-breathing and bond-stretch modes. Furthermore, the torsional dependence of the excitonic coupling was included in the latter reference.

Likewise, for purposes of comparison, standard diabatic Ehrenfest calculations were also performed. In the standard Ehrenfest method, there is only one initial electronic configuration corresponding to the left-LE expressed as a superposition of diabatic states. Each trajectory (only differing because of different initial nuclear configurations) contributes to the calculated populations of all electronic states having non-zero coupling to the initial state. Thus, Ehrenfest calculations converge with fewer trajectories than alternative methods and only 2000 trajectories were used to generate the Ehrenfest results reported here (despite being run in the diabatic representation of all 20 electronic states). The single-trajectory nature of the Ehrenfest method (as far as the electronic DOF) is, however, also the reason why it is oftentimes not particularly quantitative, as exemplified below.

The SQC/MM approach, on the other hand, includes a wide ensemble of different initial electronic configurations due to the inclusion of electronic ZPE (the Ehrenfest trajectory is the single limiting case with no ZPE within the SQC ensemble), and this variability is what allows for the statistical treatment of final electronic state outcomes which properly associates different trajectories with different final electronic states. Reflective of this, however, is that—although the SQC approach does not face the inherent exponential scaling costs of full quantum dynamics simulations—it is still more costly than a standard Ehrenfest calculation, and 30 000 trajectories were used to generate the converged SQC/MM results reported here.

Sampling of initial conditions for the simulations (ML-MCTDH, Ehrenfest, and SQC/MM) is detailed above in Sec. II. As stated, the ML-MCTDH calculations were initialized via preliminary imaginary-time propagation in the electronic subspace, and the initial electronic amplitudes obtained from this were also used to initialize the electronic DOF in the Ehrenfest calculations, whereas the electronic DOF in the SQC/MM calculations were initialized via the procedure described above involving Eqs. (20)–(22). As for the nuclear DOF, as detailed in Sec. II, the single torsion treated dynamically is initialized to θ_{10 11} = 110°; however, in the Ehrenfest and SQC/MM calculations, the ring-breathing and bond-stretch modes were initialized by two different methods: (i) the coordinates and momentum were sampled from the Wigner distribution corresponding to the ground state monomer potentials or (ii) the coordinates were initialized at exactly the minimum energy of the ground state monomer potentials with conjugate momenta taken to be zero. Both choices are presented because both have some justification: On the one hand, use of the Wigner distribution puts ZPE into the nuclear DOF which, of course, is present quantum mechanically. On the other hand, ZPE remains confined quantum mechanically, whereas classically it can and will flow amongst the various nuclear DOF without restraint. Classically, the minimum energy configuration without any destabilizing extra energy is the configuration with no ZPE.

Figure 2 shows the Ehrenfest and SQC/MM results for both initialization methods (i) and (ii) in comparison with full quantum ML-MCTDH calculations. SQC/MM results using choice (i)—the Wigner distribution (to initialize the nuclear DOF)—are seen to agree better with the ML-MCTDH results than the corresponding Ehrenfest results. In the benchmark ML-MCTDH simulation, the initially partially delocalized exciton is self-trapped on the “left” segment due to its coupling with the high frequency vibrational modes during the first 300 femtoseconds, after which the exciton hops to the center of the oligomer chain driven by planarization of the central bond torsional angle¹⁷ where it remains localized for the remainder of the simulation. In the SQC/MM simulations, all of these events are well-reproduced with reasonable accuracy. By contrast, the Ehrenfest simulations give results having the exciton more delocalized across the oligomer chain subsequent to central torsion planarization.

As for the adiabatic population dynamics shown in Fig. 2, the SQC/MM results show almost all of the initial S2 state population (∼50%) quickly transferring to the S1 state within 50 femtoseconds, and the S1 population remains dominant (>97%) afterwards and until the end of the simulation, while the S2 population remains small (<3%) and the S3 and S4 populations remain almost negligible (<1%). By contrast, the Ehrenfest approach predicts a fast but incomplete S2 → S1 transition within 50 femtoseconds, followed by a gradual decrease in the S1 population and an increase in populations of all other adiabatic states. At the end of the simulation, the S1 population is less than ∼60% and S2, S3, and S4 states have roughly the same populations. It is therefore evident that the adiabatic populations predicted by the SQC/MM simulation are in better agreement with the benchmark quantum ML-MCTDH results than the Ehrenfest simulation. The Ehrenfest method’s poor performance with regard to the adiabatic population dynamics here is likely a consequence of its well-known failure to satisfy the principle of detailed balance. The better agreement of the SQC/MM simulations with the quantum calculations suggests that it does provide a much more reasonable treatment of detailed balance, as described in Ref. 32.

Using method (ii) for initializing the nuclear DOF—minimum energy configuration with zero momentum—the agreement between SQC/MM [Fig. 2(c)] and ML-MCTDH [Fig. 2(e)] is significantly improved, in terms of torsional dynamics as well as diabatic populations where the SQC/MM results exhibit strong exciton localization throughout the duration of the simulation, as a consequence of the more restrictive sampling method. By contrast, sampling method (ii) causes the Ehrenfest results [Fig. 2(d)] to deviate more significantly from the ML-MCTDH calculations, with the exciton transferring back and forth many times across multiple monomer units, which is qualitatively inconsistent with the ML-MCTDH calculations. The torsional and adiabatic population dynamics computed with the Ehrenfest method also show qualitative differences versus the ML-MCTDH simulation.

It is interesting to speculate as to why the SQC/MM approach gives better results with sampling method (ii) than with method (i). As stated above, method (i) includes the correct quantum ZPE, but, in a classical simulation, there is no restriction on the flow of ZPE among the various vibrational modes. The results in Fig. 2 thus suggest that this unrestricted flow of ZPE destabilizes the exciton and that it is simply more consistent with a classical dynamics model not to include it. Another possibility is to include a fractional ZPE (as done with the electronic DOF in the SQC/MM approach), but this option was not explored. In any event, one also speculates that the differences in results between sampling choices (i) and (ii) would be less pronounced at higher temperatures (here the temperature is zero), where the inconsistencies between quantum and classical mechanics are generally less significant. The zero temperature case was only treated in this work to facilitate the ML-MCTDH calculations and SQC/MM benchmarking; nonzero temperature regimes are obviously more relevant to real applications and will thus be considered in future work.

Finally, note that the adiabatic population dynamics shown in Fig. 2 validates the use of only 4 adiabatic states. In particular, Fig. 2 shows that the ultrafast S2 → S1 transition within the first 50 femtoseconds is the dominant nonadiabatic transition process and that the S3 and S4 populations are nearly totally negligible during the duration of the ML-MCTDH and SQC/MM simulations. Clearly, treating just 4 adiabatic states is more than sufficient to form a complete quantitative picture of the nonadiabatic exciton dynamics here.

The SQC/MM and ML-MCTDH simulations are also consistent in illustrating that the exciton population dynamics is strongly coupled to the motion of the nuclear DOF. The time evolution of the average values of bond-stretch and ring-breathing modes at each site, calculated by averaging over the ensemble of trajectories that are windowed in the adiabatic representation at any given time t > 0, is plotted in Fig. 3—with Fig. 3(a) showing the SQC/MM simulations and Fig. 3(c) showing the ML-MCTDH results. It is evident in both calculations that the ring-breathing and the bond-stretch modes closely follow the transfer of the exciton population. When the exciton is localized in the “left” segment (monomers 1-10), the vibrational modes on the same segment fluctuate around their excited state minima in the 1D PES for each individual mode, while those on the “right” segment (monomers 11-20) fluctuate around their ground state minima. Likewise, after the exciton localizes in the center of the OT chain, the vibrational modes of the monomers near the center fluctuate around the minima of their excited state potentials, while those away from the center fluctuate around their ground state minima. This clearly indicates that electronic-nuclear coupling causes the exciton self-trapping effect, as described in Ref. 17, and it is correctly reproduced by the classical SQC/MM simulation results of Fig. 3(a) in comparison with the benchmark ML-MCTDH results of Fig. 3(c).

This analysis is confirmed by the results shown in Fig. 3(b) (SQC/MM) and Fig. 3(d) (ML-MCTDH) which differ from those shown in Figs. 3(a) and 3(c) in that the ring-breathing modes have been held fixed during the simulations (the bond-stretch modes remaining free); in these calculations, it is seen that the exciton is more delocalized (i.e., wider) versus the results in Figs. 3(a) and 3(c) where the ring breathing modes are allowed to move freely, and thus, one concludes that the self-trapping of the exciton is weaker in the absence of the ring-breathing modes. See Ref. 17 for further details regarding this effect.

Interestingly, the dynamics of the bond-stretch modes are also affected by fixing the ring-breathing modes. It is evident that in the absence of the ring-breathing modes, the bond-stretch modes take more positive values than when the ring-breathing modes are allowed to move. This reflects the intricate nature of the electronic-nuclear coupling present in the OT system. More specifically, fixing the ring-breathing modes reduces the self-trapping of the exciton population and makes the exciton more delocalized. Because the force felt by each bond-stretch mode is a weighted average of its excited and ground state PES gradients, the change in the exciton localization modifies the equilibrium geometry of the bond-stretch modes, giving rise to the more positive values. It is noteworthy that the totally classical SQC/MM approach is able to correctly reproduce these subtle effects—and the causal relationship with the dynamics of the ring-breathing modes—quite consistently with the predictions of the full quantum ML-MCTDH calculations.

The SQC/MM approach has been applied to ultrafast exciton migration dynamics in a Frenkel/site-exciton model¹⁷ of the organic semiconducting polymer oligothiophene (OT). The classical dynamics were carried out in the adiabatic representation—using the new “kinematic momentum” equations of motion (EOM) developed in Ref. 38—from which site-local (diabatic) electronic populations were extracted using a new SQC procedure developed (and validated) in this work. This new SQC procedure allows the classical simulation of exciton dynamics in Frenkel/site-exciton systems using only a few important adiabatic states, greatly enhancing the efficiency of the approach.

In this OT 20-mer system, compared with benchmark full quantum ML-MCTDH results also presented herein, the SQC/MM approach is shown to accurately describe the self-trapping of an initially partially delocalized exciton and its ultrafast migration driven by the planarization of a central torsion defect in the OT chain. It is found that the SQC/MM approach provides a more accurate description of the exciton self-trapping effect and the adiabatic populations dynamics than the traditional Ehrenfest approach, and it is speculated that this is primarily because the SQC/MM approach provides a better treatment of detailed balance versus a traditional Ehrenfest calculation. In view of these findings and methodological developments, one speculates that it would be interesting (and useful) to apply the SQC/MM approach to ultrafast exciton dynamics in semiconducting polymer systems of much larger size, possibly consisting of hundreds of monomers units, which are not easily accessible by full quantum dynamics methods. SQC/MM calculations could also be used to investigate and assess temperature effects on exciton migration in OT and similar semiconducting polymeric materials.

See supplementary material for the specification of the ground and excited state monomer potentials in Eqs. (3)–(5) and the details of the kinetic energy term in Eq. (1) as well as a description of the ML-MCTDH setup.

This work was supported by the National Science Foundation under Grant No. CHE-1464647, by the Director, Office of Science, Office of Basic Energy Sciences, Chemical Sciences, Geosciences, and Biosciences Division, U.S. Department of Energy under Contract No. DE-AC02-05CH11231, and by the Deutsche Forschungsgemeinschaft (DFG) under Grant No. BU-1032-2.

In addition, this research utilized computation resources provided by the National Energy Research Scientific Computing Center (NERSC), which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231.

The nuclear DOF in the adiabatic MM Hamiltonian [Eq. (9)] appear in terms of the usual canonical coordinates and momenta (R, P), and in terms of these variables, it is clear that when Hamilton’s equations are applied to Eq. (9), the resulting canonical EOM contain terms involving $∂ Δ P ∂ R$ and hence $D k k ′ ≡ ϕ k ∂ 2 ϕ k ′ ∂ R 2$⁠. As discussed, this can potentially be a significant drawback if the second-derivative coupling vectors, D_kk′, are not available, as is usually the case.

Fortunately, however, as shown in Ref. 38, this problem is totally (and simply) eliminated by re-expressing the EOM in terms of the “kinematic momentum” defined by

In terms of P_kin, the complete set of EOM (nuclear + electronic) corresponding to the adiabatic MM Hamiltonian [Eq. (9)]—and those used exclusively in this work—are given by

which are seen to be free of any explicit reference to second-derivative coupling vectors [particularly, Eq. (A2d)]. Notably, however, even though the second-derivative couplings do not explicitly appear, there is no approximation involved; it turns out, as described in Ref. 38, that the effects of the second-derivative couplings are still fully accounted for.

Lastly, regarding the effective potential term V_eff appearing in Eq. (A2d): Applying Hamilton’s equations to Eq. (9), this term would be given by

In practice, however, a symmetrized version of the MM Hamiltonian is almost always employed, to which the EOM of Eq. (A2) still correspond so long as the effective potential term is taken to be

Equations (A2) and (A4) thus constitute the operative EOM for all calculations here. (Again, for the details, see Ref. 38.)

In order to validate the SQC procedure set forth in Sec. II B 3 for extracting diabatic electronic populations from adiabatic MM dynamics calculations, a similar calculation is shown here for the set of spin-boson model problems previously treated in Ref. 30 by the standard diabatic SQC/MM approach. These results, as shown in Fig. 4, are plotted against the standard SQC/MM results. Specifically, to demonstrate the diabatic population extraction procedure, the new results were generated by initializing trajectories in the diabatic representation, converting them to the adiabatic representation, running the adiabatic dynamics to time t, windowing the adiabatic electronic variables at time t, including the off-diagonal density matrix element, and transforming the results via Eq. (23) to compute the diabatic populations. Figure 4 verifies that this new procedure does indeed do a reasonable job of reproducing SQC/MM results calculated directly in the diabatic representation (which were shown in Ref. 30 to be in good agreement with the benchmark QM results).