In the last several years, a symmetrical quasiclassical (SQC) windowing model applied to the classical MeyerMiller (MM) vibronic Hamiltonian has been shown to be a simple, efficient, general, and quiteaccurate method for treating electronically nonadiabatic processes at the totally classical level. Here, the SQC/MM methodology is applied to ultrafast exciton dynamics in a Frenkel/siteexciton model of oligothiophene (OT) as a model of organic semiconductor polymers. In order to keep the electronic representation as compact and efficient as possible, the adiabatic version of the MM Hamiltonian was employed, with dynamical calculations carried out in the recently developed “kinematic momentum” representation, from which site/monomerspecific (diabatic) excitation probabilities were extracted using a new procedure developed in this work. The SQC/MM simulation results are seen to describe coherent exciton transport driven by planarization of a central torsion defect in the OT oligomer as well as to capture exciton selftrapping effects in good agreement with benchmark quantum calculations using the multilayer multiconfiguration timedependent Hartree approach. The SQC/MM calculations are also seen to significantly outperform the standard Ehrenfest approach, which shows serious discrepancies. These results are encouraging, not only because they illustrate a significant further application of the SQC/MM approach and its utility, but because they strongly suggest that classical mechanical simulations (with the potential for linear scaling efficiency) can be used to capture, quantitatively, important dynamical features of electronic excitation energy transfer in semiconducting polymers.
I. INTRODUCTION
Exciton dynamics, or excitation energy transfer (EET), in conjugated semiconducting polymers plays an important role in determining the efficiency of organic electronic materials.^{1} 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 210 monomer units.^{3} 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.^{13} Additionally, the formation and relaxation of local exciton ground states (LEGS),^{14–16} which are subject to exciton selftrapping,^{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 highdimensional systems exhibiting strong electronphonon coupling. Classical investigations of the Ehrenfest type^{15} and surface hopping variety^{19} have been undertaken. They are, however, of limited accuracy as far as the treatment of electronphonon correlations is concerned. Full quantum dynamical studies using efficient multiconfigurational methods, i.e., the multiconfiguration timedependent Hartree (MCTDH) method^{20,21} and its multilayer (MLMCTDH) variant,^{22–24} are accurate but have generally been restricted with regard to the system size.^{25–28} For systems composed of up to 2030 sites, recent MLMCTDH simulations support experimental observations that the relaxation of geometric defects such as torsional modes drives the exciton transfer across monomer units^{17} 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 modestsized 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 wellknown 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^{37} to include the treatment of electronically nonadiabatic processes—i.e., those molecular processes which involve multiple electronic states/BornOppenheimer 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^{29} 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^{30} to the electronic DOF of the MM Hamiltonian,^{39} in combination with a corresponding matching of the zeropoint 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^{30} and photochemistry models, condensedphase spinbosontype models of dissipative electronic dynamics,^{30} electron transfer reactions,^{31} Frenkeltype siteexciton models of lightharvesting complexes,^{34} etc. In all these application contexts, the SQC/MM methodology is seen capable of quantitatively describing electronic quantum coherence, from regimes of extremely strong^{37} to extremely weak^{36} coupling between the electronic states, and including the decoherence 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^{35} 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^{36} which, though developed to deal with the weakcoupling limit, is found superior to the original square histogram SQC model in all coupling regimes; and, most recently, (iii) a new “kinematic momentum” representation^{38} for adiabatic MM vibronic dynamics which eliminates—analytically, without approximation—the need to calculate secondderivative couplings (as in the adiabatic multistate/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 BornOppenheimer potential energy surfaces (PESs) are generated “onthefly” 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(3hexylthiophene) (P3HT), a frequently employed donor material in organic photovoltaics.^{1} Current theoretical investigations of exciton dynamics in such species often rely on the use of model Hamiltonians, and in recent MCTDH and MLMCTDH studies,^{26–28} a Frenkel siteexciton model Hamiltonian for OT type systems was derived from electronic structure calculations for suitable oligomer species,^{40} 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 lowlying adiabatic electronic states, the remaining states not requiring explicit treatment.
For benchmarking, SQC/MM calculations are compared with exact quantum calculations obtained via the multilayer multiconfiguration timedependent Hartree (MLMCTDH) 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.^{15} The full quantum benchmark results (obtained via MLMCTDH) are closely analogous to the simulations of Ref. 17, which provide evidence for the dynamics of selftrapping and torsionally induced exciton migration in a lowtemperature regime. Below, it is shown that the SQC/MM simulation results agree well with the benchmark MLMCTDH calculations and significantly outperform the Ehrenfest simulations.
The remainder of the paper is organized as follows: In Sec. II, the model Frenkel/siteexciton OT Hamiltonian is presented, as well as the details of the MLMCTDH, 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.
II. THEORY
The exciton dynamics study presented in this paper centers on the treatment—by both quantum and classical mechanics—of a Frenkel/siteexciton 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 MLMCTDH^{22–24} methodology, and relevant details of these calculations will be discussed as they arise.
A. Oligothiophene model
1. The Frenkel siteexciton Hamiltonian
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 singlemonomer excitations in the directproduct 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/sitedependent potentials $ { V k site } $, each corresponding to an electronic configuration localizing an excitation on the kth of the F monomers, a stateindependent bath component V_{bath}, and excitonic (diabatic) couplings $ { V k \u2009 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 sitelocal ringbreathing modes, X ≡ {X_{k}}, along with two types of “correlated” modes that connect neighboring monomer sites: interring bondstretch modes, Y ≡ {Y_{k}_{k+1}}, and torsional modes, θ ≡ {θ_{k}_{k+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 ringbreathing mode (X_{k}), two bondstretches (Y_{kk−1}, Y_{kk+1}), and torsions (θ_{kk−1}, θ_{kk+1}) shared with adjacent monomers.
The state/sitedependent potentials in Eq. (2) are expressed in terms of the given nuclear modes (X, Y, θ) as
with the groundstate 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 pointwise 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 highfrequency ringbreathing and bondstretch 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 MLMCTDH calculations to a reasonable computational effort, only the intermonomer torsion at the middle of the OT segment, θ_{10 11}, is treated dynamically (which, as explained below, serves to break the segment into two 10mer “spectroscopic units”). Accordingly, only θ_{10 11} is coupled to the harmonic bath, which is thus given by
where, as usual (and assuming massweighted coordinates), the bath couplings {c_{j}} are adapted to an Ohmic spectral density, $ c j = \omega j ( 2 / \pi ) \u2009 \eta \u2009 \Delta \omega $ 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.
2. Initial conditions
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 selftrapping 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 unitdefining 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 nonplanar configuration, breaking the πconjugation at θ_{10 11}, effectively defining two spectroscopic units as left and right 10mer segments. Also, as stated above, θ_{10 11} is the only torsional bond angle treated dynamically; all others are locked at their groundstate equilibrium geometries of θ_{eq} = 160°.
Within the singleexcitation Frenkel manifold, the lowestlying 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 (leftLE) and this, of course, must be properly set up for the dynamical simulations.
For the full quantum MLMCTDH calculations, the initial leftLE electronic configuration is established via imaginarytime propagation in the electronic subspace, prior to the realtime quantum propagation.^{17,26} This corresponds to the lowest adiabatic electronic eigenstate of the left 10mer OT segment.
The electronic amplitudes emerging from the imaginarytime 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 groundstate potential V_{0} 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 groundstate equilibrium geometry of θ_{eq} = 160°. Thus, the initial nuclear conditions correspond to a vertical (FrankCondon) excitation into an electronic state with delocalization broken by the torsional defect.
B. The SQC/MM approach
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 MeyerMiller (MM) Hamiltonian and (ii) the symmetrical quasiclassical (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.
1. The MeyerMiller (MM) Hamiltonian
The MeyerMiller (MM) Hamiltonian^{39} 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 offdiagonal 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 zeropoint energy (ZPE) parameter, originally taken to be $ 1 2 $ (the quantum value) but, in the SQC approach, always taken to be $\u2248 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^{43} that if this Hamiltonian was treated quantum mechanically—i.e., as an operator in the timedependent 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 welldefined classicallimit of the true multistate vibronic quantum dynamics.
The timedependent Schrödinger equation for multiple electronic states may also of course be written in an adiabatic representation—i.e., in terms of the BornOppenheimer PESs {E_{k}(R)} and first and secondderivative 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 offdiagonal 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 lowestlying 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 secondderivative 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 onthefly ab initio dynamics calculation, as secondderivative 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.
2. Symmetrical quasiclassical (SQC) quantization
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 quasiclassical (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 prefactors 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 actionangle oscillator variables {n_{k}, q_{k}} (in terms of which the windowing/quantization is done) is
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 timedependent 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 actionangle 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 offdiagonal 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 phaseangles 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 \u2032 ) = ( 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_{0}, R_{0}) for the G nuclear DOF amounts to Monte Carlo evaluation of the phasespace average
for all k, k′ ∈ [1, F], followed by renormalization via
In other words, in an SQC/MM simulation, trajectories are

initialized via importance sampling: ρ (P_{0}, R_{0}) is used to sample the G nuclear DOF, and the F electronic actions are sampled according to W_{i}(n_{0}) with initial angles q_{0} chosen randomly in the interval [0, 2π] and (n_{0}, q_{0}) converted to Cartesian (p_{0}, x_{0}) according to Eq. (12),

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

collected at time t via W_{kk′}(n, q) [Eq. (17)]—after conversion to actionangle variables (n_{t}, q_{t}) according to Eq. (11)—to calculate $ { \rho \u0303 k k \u2032 \u2190 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 $ { \rho \u0303 k k \u2032 \u2190 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 offdiagonal density matrix elements.
3. Boundary conditions for the SQC/MM simulations
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 24 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 sitelocal 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 leftLE 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 (Frenkelexciton) diabatic representation describes this sitespecific localization of electronic excitation.
a. Setting up the initial leftLE electronic configuration.
As described in Sec. II A, the 20mer OT system is to be initialized with a torsional defect at its center, splitting it into a pair of 10mer “spectroscopic units,” whereby the initial electronic excitation is the (partially) localized exciton ground state of the left segment—i.e., the initial leftLE. 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 ringbreathing and bondstretch 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 MLMCTDH 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 leftLE configuration can be established for the electronic DOF. Because the dynamics is to be run in the adiabatic representation, the protocol requires the nuclear coordinatedependent adiabatictodiabatic transformation matrix obtained from diagonalizing the electronic Hamiltonian Ĥ_{elec}(R) [Eq. (2)]. This unitary transformation is designated
where $ { \varphi 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 braket refers to the underlying diabatic/Frenkel basis.
The leftLE initial condition, however, is not the ground adiabatic electronic state of the full Ĥ_{elec} but rather the lowestlying electronic state with excitation localized on the left segment—i.e., it is the ground adiabatic state of the left 10mer, which roughly corresponds to a coherent superposition of the two lowestlying adiabatic states of the full 20mer. This electronic state may be determined by blockdiagonalizing Ĥ_{elec} in the diabatic basis with the coupling(s) between left and rightsegments set to zero,
The corresponding coordinatedependent unitary matrix—transforming from this coordinatedependent LE electronic basis $ { \varphi k L E } $ [the eigenvectors of Eq. (20)] to the diabatic/Frenkel basis—is designated as
The portion of the lowestlying 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 leftLE electronic variables to the adiabatic representation for dynamics propagation. Having sampled the initial nuclear configuration from ρ (P_{0}, R_{0}) and determined the two transformation matrices of Eqs. (19) and (21), the initial sampled electronic variables (n_{0}, q_{0}) are converted to Cartesian (p_{0}, x_{0}) via Eq. (12) and then further transformed to the adiabatic basis by sequentially applying the foregoing unitary transforms,
the first transforming from leftLE 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 \u2009 p k ) $ [Eq. (13)].
This completes the selection of the initial trajectory variables [step 1 in evaluating Eq. (18)] consistent with the leftLE 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 leftLE 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.
b. Extracting the final diabatic electronic populations from adiabatic MM dynamics.
The final step of the SQC protocol—after timeevolving the electronic and nuclear DOF—is to window the final electronic actionangle 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 20mer OT chain and, of course, this is naturally tied to the diabatic/Frenkel sitelocal electronic basis. Accordingly, some protocol must be devised for extracting sitelocal 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 leftLE. 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 adiabatictodiabatic unitary matrix of Eq. (19)—is coordinatedependent, 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 adiabatictodiabatic transformation matrix of Eq. (19). That is, each trajectory propagated in the adiabatic basis makes a contribution to the diabatic/sitelocal 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 timeevolved actions and angles (n_{t}, q_{t}) on the RHS are adiabatic, but the LHS is the desired diabatic/Frenkelexciton 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,^{44} and these methods also rely on computing offdiagonal 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 spinboson problems, where it is indeed shown that trajectory propagation and windowing in the adiabatic electronic basis, followed by extraction of diabatic populations via a coordinatedependent adiabatictodiabatic 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 leftLE 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 leftLE 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.
III. RESULTS
A. Benchmarking versus full quantum and classical Ehrenfest calculations
Benchmark quantum dynamical calculations were performed using the MLMCTDH method^{23,24} (Heidelberg MCTDH package^{45}) in the diabatic representation of Eqs. (1) and (2) and comprising a sevenlayer 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 ringbreathing and bondstretch 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 leftLE 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 nonzero 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 singletrajectory 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 (MLMCTDH, Ehrenfest, and SQC/MM) is detailed above in Sec. II. As stated, the MLMCTDH calculations were initialized via preliminary imaginarytime 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 ringbreathing and bondstretch 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 MLMCTDH calculations. SQC/MM results using choice (i)—the Wigner distribution (to initialize the nuclear DOF)—are seen to agree better with the MLMCTDH results than the corresponding Ehrenfest results. In the benchmark MLMCTDH simulation, the initially partially delocalized exciton is selftrapped 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^{17} where it remains localized for the remainder of the simulation. In the SQC/MM simulations, all of these events are wellreproduced 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 MLMCTDH 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 wellknown 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 MLMCTDH [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 MLMCTDH calculations, with the exciton transferring back and forth many times across multiple monomer units, which is qualitatively inconsistent with the MLMCTDH calculations. The torsional and adiabatic population dynamics computed with the Ehrenfest method also show qualitative differences versus the MLMCTDH 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 MLMCTDH 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 MLMCTDH 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.
B. SQC/MM demonstration of exciton selftrapping
The SQC/MM and MLMCTDH 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 bondstretch and ringbreathing 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 MLMCTDH results. It is evident in both calculations that the ringbreathing and the bondstretch modes closely follow the transfer of the exciton population. When the exciton is localized in the “left” segment (monomers 110), 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 1120) 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 electronicnuclear coupling causes the exciton selftrapping 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 MLMCTDH results of Fig. 3(c).
This analysis is confirmed by the results shown in Fig. 3(b) (SQC/MM) and Fig. 3(d) (MLMCTDH) which differ from those shown in Figs. 3(a) and 3(c) in that the ringbreathing modes have been held fixed during the simulations (the bondstretch 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 selftrapping of the exciton is weaker in the absence of the ringbreathing modes. See Ref. 17 for further details regarding this effect.
Interestingly, the dynamics of the bondstretch modes are also affected by fixing the ringbreathing modes. It is evident that in the absence of the ringbreathing modes, the bondstretch modes take more positive values than when the ringbreathing modes are allowed to move. This reflects the intricate nature of the electronicnuclear coupling present in the OT system. More specifically, fixing the ringbreathing modes reduces the selftrapping of the exciton population and makes the exciton more delocalized. Because the force felt by each bondstretch 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 bondstretch 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 ringbreathing modes—quite consistently with the predictions of the full quantum MLMCTDH calculations.
IV. CONCLUSION
The SQC/MM approach has been applied to ultrafast exciton migration dynamics in a Frenkel/siteexciton model^{17} 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 sitelocal (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/siteexciton systems using only a few important adiabatic states, greatly enhancing the efficiency of the approach.
In this OT 20mer system, compared with benchmark full quantum MLMCTDH results also presented herein, the SQC/MM approach is shown to accurately describe the selftrapping 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 selftrapping 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.
SUPPLEMENTARY MATERIAL
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 MLMCTDH setup.
ACKNOWLEDGMENTS
This work was supported by the National Science Foundation under Grant No. CHE1464647, 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. DEAC0205CH11231, and by the Deutsche Forschungsgemeinschaft (DFG) under Grant No. BU10322.
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. DEAC0205CH11231.
APPENDIX A: “KINEMATIC” ADIABATIC REPRESENTATION AND EQUATIONS OF MOTION (EOM)
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 $ \u2202 \Delta P \u2202 R $ and hence $ D k k \u2032 \u2261 \varphi k \u2202 2 \varphi k \u2032 \u2202 R 2 $. As discussed, this can potentially be a significant drawback if the secondderivative 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 reexpressing 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 secondderivative coupling vectors [particularly, Eq. (A2d)]. Notably, however, even though the secondderivative couplings do not explicitly appear, there is no approximation involved; it turns out, as described in Ref. 38, that the effects of the secondderivative 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
APPENDIX B: VALIDATION OF THE SQC DIABATIC POPULATION EXTRACTION PROCEDURE OF SEC. II B 3
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 spinboson 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 offdiagonal 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).