# Quantum Simulations of Charge Separation at a Model Donor-Acceptor Interface: Role of Delocalization and Local Packing

Authors: Allen Kelley; Kush Patel; Eric R. Bittner
Journal: Advances in Condensed Matter Physics
 (2018)

---

## Abstract

Organic Polymer-based photovoltaic systems offer a viable alternative to more standard solid-state devices for light-harvesting applications. In this study, we investigate the electronic dynamics of a model organic photovoltaic (OPV) heterojunction consisting of polyphenylene vinylene (PPV) oligomers and a[6,6]-phenyl C61-butyric acid methyl ester (PCBM) blend. Our approach treats the classical molecular dynamics of the atoms within an Ehrenfest mean-field treatment of the π-π⁎ singly excited states spanning a subset of donor and acceptor molecules near the phase boundary of the blend. Our results indicate that interfacial electronic states are modulated by C=C bond stretching motions and that such motions induce avoided crossings between nearby excited states thereby facilitating transitions from localized excitonic configurations to delocalized charge-separated configurations on an ultrafast time-scale. The lowest few excited states of the model interface rapidly mix allowing low frequency C-C out-of-plane torsions to modulate the potential energy surface such that the system can sample both intermolecular charge-transfer and charge-separated electronic configurations on sub-100 fs time scales. Our simulations support an emerging picture of carrier generation in OPV systems in which interfacial electronic states can rapidly decay into charge-separated and current producing states via coupling to vibronic degrees of freedom.


---

## Body

## 1. Introduction

Advances in both materials and device fabrication have led to the development of highly efficient organic polymer-based photovoltaic cell (OPV) in which the power conversion efficiency is in excess of 10-11% under standard solar illumination and efficiencies as high as 12% in multijunction OPVs [1]. This increase in power conversion efficiency indicates that mobile charge carriers can be efficiently generated and collected in well-optimized devices; however, the underlying photo-physical mechanism for converting highly bound molecular (Frenkel) excitons into mobile and asymptotically free photocarriers remains elusive in spite of vigorous, multidisciplinary research activity [2–18].It is generally recognized that the local microstructure of the interface between the donor and acceptor materials composing a bulk heterojunction cell plays a central role in determining the overall efficiency of an OPV device. Generally speaking, charge separation in OPV materials occurs when excitons near the donor/acceptor (DA) interface dissociate into either charge-transfer or charge-separated states.(1)D∗:A⟶D+:A-orD+⋯A-These two states are distinguished by whether or not the positively charged (hole) on an electron donor molecule is directly adjacent to the acceptor anion, termed a charge-transfer state (CT), or if the two charges are separated by one or more molecular units, termed a charge-separated state (CS). In the former (CT) case, the two charges are pinned to the D : A interface by Coulombic forces which are only weakly screened by the low dielectric nature of conjugated organic systems. Curiously, recent ultrafast spectroscopic measurements on optimized polymer : fullerene-based OPV systems indicate that charge carriers begin to appear ≤ 100-fs following photoexcitation [2, 12, 15, 19–22]. Moreover experiments by Gélinas et al., in which Stark-effect signatures in transient absorption spectra were analysed to probe the local electric field as charge separation proceeds, indicate that electrons and holes separate by as much as 40 Å over the first 100 fs and evolve further on a picosecond time scales to produce unbound and freely mobile charge pairs [4].It is generally recognized the electronic states at DA interfaces are delocalized over multiple molecular units and may well extend deeply into the bulk. As a result, the spatial extent of these states presents a profound challenge to even state-of-the-art ab initio quantum chemical methods. Moreover, one can not ignore the role of structural fluctuations and intermolecular vibrations in modulating the electronic energies. Recently, Bittner and Silva proposed that energy level fluctuations due to thermal motions can lead to strong quantum mechanical coupling between excitonic and charge-separated states [23]. Central to this theory is the fact that the golden-rule transition rates between states can be determined by examining energy gap fluctuations. This model unifies a number of seeming disparate observations but lacks a desired atomistic basis.Here we employ a hybrid molecular mechanics/quantum mechanics (MM/QM) approach to sample the electronic excited states of a model bulk heterojunction system. Our approach provides an explicit, semiempirical description of the excitedπ-electronic states of select molecular units near the interface between phase segregated donor and acceptor molecular units. Our simulations indicate that interfacial electronic states are modulated by C=C bond stretching motions and that such motions induce avoided crossings between nearby excited states thereby facilitating transitions from localized excitonic configurations to delocalized charge-separated configurations on an ultrafast time-scale.


## 2. Results

Our results build upon a model of charge separation in donor : acceptor blends that we recently reported upon [24]. In this study, we considered a model blend of poly(p-phenylene vinylene) (PPV) donor oligomers and phenyl-C61-butyric acid methyl ester (PCBM) and defined the “π-active” region by selecting (by hand) a single PCBM and three nearby PPV chains to assess the penetration of a localized electronic state into the bulk region. In this paper we model a phase separated interface between PCBM and PPV oligomers. While the semiempirical treatment does facilitate including up to a few hundred π electrons, we do need to be judicious in choosing the “π-active” region.Figure1 shows representative molecular configurations of the three cases studied in this paper. Each corresponds to a periodic simulation cell in the xy plane following equilibration at 100 K and 1 atm pressure and they are representative of typical interfacial configurations. In each figure, the red and blue colored spheres represent atoms included in the quantum chemical description. In Case  A, we selected 2 interfacial PCBMs and 3 nearby π-active PPV oligomers that penetrate into the bulk polymer region, including a total of 230 C atoms. In Case B, we selected 3 PCBM and 3 nearby PPV oligomers expanding the π-active molecules that form the interfacial heterojunction, including a total of 288 C atoms. Cases  A and B differ solely in the initial selection of the PPV oligomers included in the quantum mechanical subspace. For comparison, in Case  C we selected 1 PCBM and 3 PPV oligomers that penetrate into the bulk. In all cases, we employ periodic boundary conditions in xyz as in our previous study [24]. This simulation was set up to be very similar to case A, only adding a single PCBM molecule to the interfacial region.Figure 1
Snapshot of the three MD simulation cells: A, B, and C. The red and blue highlighted molecules denote theπ-active units in our system.
(a) (b) (c)These three cases are certainly not an exhaustive study of all possible ways that an excitation can localize at an interface and in general one should include all the C atoms as being part of theπ system as to include the possibility for both localized and delocalized excitations. However, such study is not computationally feasible and we select these to test the extent that low-lying excitations can penetrate into the bulk region, how this penetration may facilitate charge separation, and the extent that this plays into the mixing between charge-transfer, excitonic, and charge-separated configurations.Over the course of eight 10 ps simulations, the lowest lying excitation samples a variety of electronic states ranging from localized excitons to charge-separated, charge-transfer, and delocalized states. Figure2 shows the charge density for four typical eigenstates encountered over the course of a given simulation. We characterize these as follows:(i)
Excitonic (EX) states characterized as having>50% of the electron/hole density on a single molecule. The exciton state populates ≈ 58% of the states for simulation A, ≈32% for simulation B, and ≈52% for simulation C.(ii)
Charge-Transfer states (CT) characterized as having>50% of the electron/hole density occupying adjacent molecules. The charge-transfer state populates ≈ 14% of the states for simulation A, ≈19% for simulation B, and ≈12% for simulation C.(iii)
Charge-Separated (CS) states characterized as having>50% of the electron/hole density occupying a PCBM and a PPV separated by a single molecule. The charge-separated state are populated ≈ 26% for simulation A, ≈49% for simulation B, and ≈26% for simulation C.(iv)
Delocalized states (DLOC), characterized as having<50% of the electron/hole density on a single PPV or PCBM molecule.Figure 2
The dominant adiabatic states from simulations A (c, d) and B (a, b) as shown in Figure1 The red and blue numbers denote the electron/hole density as a percent on the indicated molecule. The states correspond to the x-axis in Figure 3. The four snapshots shown represent typical states for our systems. (a, b) from left to right present an exciton located on the PCBM molecule and a charge-transfer state. (c, d) from left to right show a charge-separated state and a partly delocalized state.
(a) (b) (c) (d)Figure 3
Histogram plot of the populations of adiabatic states in systems A, B, and C taken every fs. Each cluster represents a different classification of the adiabatic states as shown in Figure2.For the most part, the lowest eigenstates in simulations A and B do fall into the categories above. In simulation C, however, the lowest electronic excited state spends about 10% of its time as an eigenstate that cannot be neatly binned into the above categories and corresponds to a fairly equal mixing between excitonic and charge-transfer states and excitonic and charge-separated states.The striking difference between the blend (C) and the phase separated cases (A and B) is that the stability of localized excitons is very much determined by the proximity between donor and acceptor molecules. In both A and C, theπ electrons spend over half their time in excitonic configurations and less than 30% of their time in charge-separated configurations. In both A and C, the 3 PPV chains included in the π system are stacked into the bulk region away from the PCBMs. In contrast in B, charge-separated configurations, defined as a buildup of negative charge on the more distant PCBM, dominate. Only rarely in all cases was the excited state considered as “delocalized.”The energy of the lowest excited state following excitation att=0 fs for simulations A and B is shown in Figure 4(a). After excitation there is very little energetic relaxation in all of the systems simulated. The simulations appear to cycle through many adiabatic states in a short period of time leaving the impression of a weak electron-phonon coupling. This can be rationalized as the electron/hole density often delocalize over multiple molecules and many conjugated C-C bonds. Another striking effect of the systems is the large number of avoided crossings that occur between the lowest lying states. There is also a 20 fs oscillation in the CI energies, driving the systems excited states into many regions of strong coupling. The 20 fs oscillation also appears in the autocorrelation and the Fourier transform of the gap energies, contributing the C=C bond stretching modes around ≈1600 cm−1. The oscillation is contributed to small thermally activated fluctuations within the simulation, showing that even at 100 K the thermal fluctuations possess sufficient energy to bring these states into regions of strong electronic coupling.Figure 4
(a) Excitation energies for simulations A and B for the first 200 fs following excitation to the lowest SCI state. Each colored line corresponds to the excitation energy from the ground electronic state to one of the lowest 5 singlet excited states. (b) Histogram distributions of the 5 lowest excitation energy levels over a 200 ps simulation. Throughout the simulation, colors of each histogram correspond to the time-histories plotted in (a).
(a) (b)In Figure4(b) we show a histogram of the 5 lowest excitation energies accumulated over 40 ps of simulation time following promotion to the lowest excited singlet state. The energies take a Gaussian distribution around their respective means and the distributions overlap. On average there is an energydifference between 50 and 70 meV between the first and second excited state energies. This suggests that, on average, there is a dense nearly continuous spectrum of charge-separated excited states that can easily be brought into strong electronic coupling by small fluctuations the system brought about by nuclear motion.In other words, even though the lowest excited state generally can be considered to be an interfacial charge-transfer state, small nuclear motions rapidly mix electronic configurations allowing such states to become delocalized charge-separated states.We next consider the origins of the energy fluctuations evidenced in Figure4(a). While we only show two 200 fs segments of eight 10 ps simulations over this period, one can see that the excitation energies are modulated and cover a small range. The autocorrelation plots of the band gap energies, shown in Figure 5, show that the correlation times for the three simulations are very short ≈ 8 fs meaning that the system changes rapidly enough that the oscillations observed are independent of one another. By taking the Fourier transform of the gap the IR active modes that contribute to the modulation of the excitation energies inside of the systems are found as shown in Figure 5. The modulation of the excitation energy appears to be heavily dependent upon the torsion, C=C, and C-H stretching modes. In each of the plots three distinct regions can be seen, the low frequency torsional modes occur between 200 and 500 cm−1, the C=C stretching modes occur between 1300 and 1800 cm−1, and the C-H stretching modes occur between 2800 and 3300 cm−1. We conclude that small-scale vibronic fluctuations in the molecular structures and orientations produce significant energetic overlap between different adiabatic states to drive the system from purely excitonic to purely charge-transfer on a rapid time-scale. This is evidenced in the progression of the excitation energies, as small fluctuations in these modes can easily bring the excited states into strong coupling regimes.Figure 5
(a–c) Autocorrelation of the energy gap between first and second excited states for simulations A, B, and C. The correlation times for all three simulations is 8 fs. (d–f) These plots show the active IR active modes that contribute to the modulation in the CI energies of the systems. Three distinct regions are visible in each plot. The low frequency torsional modes occur between 200 and 500 cm−1, the C=C stretching modes occur between 1300 and 1800 cm−1, and the C-H stretching modes occur between 2800 and 3300 cm−1.
(a) (b) (c) (d) (e) (f)According to the model outlined below the mean(Δ0) and variance (V-) of the energy gap distributions shown in Figure 4 can be used as input to estimate the state to state transition rate for a two-level system. We take Td-1≈V¯/ħ as an estimate of the decoherence time and we introduce τ as the natural lifetimes of each state. The results are shown in Table 1. The estimated transition rates are consistent with the observations that the systems rapidly sample a wide number of possible configurations over the course of the molecular dynamics simulation. On average, the state to state couplings of 56 meV for simulation A and 48 meV for simulation B are comparable to the average energy gaps between the lowest excited states. The strong electronic coupling allows for rapid transitions; however, larger couplings also imply shorter electronic decoherence times, effectively quenching the ability of charges to separate by tunneling.Table 1
Estimated excited state(i→f) transition rates and vibronic couplings for simulations A, B, and C.
Transition Δ o (eV) 〈 V 〉 (eV) T d (fs) k - 1 (fs) A 1→2 0.070 0.050 13.16 19.5 1→3 0.12 0.056 11.75 32.9 1→4 0.16 0.063 10.45 38.9 B 1→2 0.050 0.042 15.67 18.9 1→3 0.09 0.049 13.43 29.3 1→4 0.12 0.052 12.66 40.0 C 1→2 0.11 0.057 11.41 21.8 1→3 0.18 0.076 8.66 28.6 1→4 0.26 0.076 8.70 50.0


## 3. Discussion

We present here the results of hybrid QM/MM simulations of the excited states of model PPV/PCBM heterojunction interfaces. Our results indicate that varying the blend ratio and placement of the molecules comprising the heterojunction greatly affect the distribution of states yet have little effect upon the rate constants of the system. We also propose that thermal noise can rapidly change the character of the lowest lying excited state from purely excitonic to charge-separated on a time-scale of sub-100 fs.Simulations A and C have a very similar placement of molecules, only differing in that simulation A adds a PCBM molecule to the heterojunction. The addition of the PCBM only slightly changes the distribution of states as seen in Figure4. The exciton states continue to be the most favoured state inside the system, even slightly increasing in probability, while the delocalized states slightly decrease in probability. Simulation B completely changes the heterojunction, placing three PPV and two PCBM molecules at the interface as seen in Figure 1. The states generated by simulation B are radically different from those seen in simulation A and C as seen in Figure 4. The probability of finding the system in the exciton state is dramatically reduced while the charge-separated state becomes predominant. This result is quite interesting as it highlights that the complexity of simulating heterojunctions resides not only in the size of the system but also on how the donor/accepter interface is chosen.All three systems start with an exciton localized on the PCBM and dissociating into a charge-transfer state with the hole (or electron) delocalized over multiple polymer units before localizing to form charge-separated states. There is a wide range of electronic states tightly clustered within a small energy band, allowing small changes in local bond lengths to have a dramatic role in modulating the electronic couplings between excited states. We speculate that the dramatic shift in population seen in simulation B can be caused by disorder in the PPV molecules reducing the band gap by 20 meV. The PPV molecules comprising the interface region undergo large distortions in the C-C torsion angles allowing the molecules to cycle through a larger range of configurations inside of a short time interval. The presence of moreπ active PPV molecules at the interface also appears to lead to more avoided crossing regions and the ability of the system to more efficiently dissociate excitons into charge-transfer and charge-separated states to a distance to where their Coulombic attraction is comparable to the thermal energy. While the finite size of our system prevents further dissociation of the charges, the results are suggestive that such interstate crossing events driven by bond-fluctuations can efficiently separate the charges. The results presented here corroborate recent ultrafast experimental evidence suggesting that free polarons can form on an ultrafast time scales (sub-100 fs) and that thermally activated low frequency torsional modes are key in effective electron hole separation in PPV/PCBM heterojunctions.


## 4. Methods

### 4.1. Simulation Methods

Our simulations employ a modified version of the TINKER molecular dynamics (MD) package [25] in which the MM3 [26] intramolecular bonding parameters are allowed to vary with the local π-electronic density as described by a Parisier-Parr-Pople (PPP) semiempirical Hamiltonian [27, 28]. Similar approaches have been described by Lobaugh and Rossky [29] and Muller and Warshel [30] to include electronic dynamics into an otherwise classical force field description.We assume that the nuclear dynamics of all atoms in the system evolve according to a classical description(2)mir¨it=Firt,Pt,where Firt,Pt is the force acting on atom i at time t and Pt is the instantaneous electronic density for the π-electrons. The Fock operator for the π electrons is given by(3)Fμμ=Wμ+12Pμμ∘μμ∣μμ+∑ρ≠μPρρ-ZρVrμρ,Fμν=βμνrμν-12Pμν∘μμ∣νν,where Wμ is the valence state potential for a given atom type and we make the zero-differential overlap (ZDO) assumption to write the electronic repulsion integrals as(4)μν∣ρσ=μμ∣ρρδμνδρσ.The Greek indices correspond to atom-centred basis orbitals. Note that we specifically distinguish P∘ as the ground state (Hartree-Fock) electronic density. Both the hopping integrals βμν and Coulomb repulsion integrals are determined using empirical formulas and depend upon the intermolecular bonding between atoms, accounting for nonplanarity. The molecular orbitals are eigenstates of F(5)Fψn=ϵnψn,where ψn is given as a linear combination of the atom-centred orbitals(6)ψn=∑μcnμχμ.From the cnμ coefficients, we obtain the ground state density matrix(7)Pμν∘=2∑nocccnμcnν.The electronic energy is given by(8)EHF=Tr⁡FP∘and is minimized by imposing the condition(9)F,P∘=0.Electronic π→π∗ excitations are subsequently introduced by using configuration interaction (singles) (CI-S) theory where we write the excited state wave function as a linear combination of elementary excitations between occupied and unoccupied molecular orbitals(10)Ψk=∑h∈occ∑e∈unoccCkehψhψe=∑μν∑ehCkehceμchνχνχμ=∑μνCkμνχνχμ,where the {Ckμν} coefficients give the electron/hole amplitude in the atom-centred basis. These are eigenstates of the CI-S Hamiltonian and we used a total of 10 occupied and 10 unoccupied Hartree-Fock molecular orbitals to construct the electron/hole configurations for the CI calculations. The excited state bond charge density matrix is constructed by assuming that electron densities are added to the unoccupied orbitals and hole densities are subtracted from the Hartree-Fock ground state; namely,(11)Pμν=Pμν∘+Pμνe-Pμνh,where(12)Pμνe=∑σCkμσCkνσPμνh=∑σCkσμCkσν.In the former, we sum over the occupied orbitals (holes) and in the latter we sum over unoccupied orbitals. Pμν is interpreted as the “bond charge” matrix whereby the diagonal elements give the net π-electron count on site μ and the off-diagonal Pμν terms give the π-bond order between sites μ and ν. Consequently, excitations systematically move electron density from atom to atom and change the local bond orders.Within the force field, bond force field parameters are modified by assuming a linear interpolation between single and double bonds. For example, we write the bond force constant between C atomsi and j as(13)k=k2-δk1-Pij,where k2 is the force constant for a C=C double bond when Pij=1, k1 is the bond force constant for a C-C single bond when Pij=0, and δk=(k2-k1). In a similar definition, we interpolate the equilibrium bond lengths between single- and double bonds as(14)l=l2+δl1-Pij.This treatment is not limited to C-C bonds but is applied to all heavy atoms contributing to the π system. While the MM3 parameter set was originally intended to account for the electronic ground state, we assume that the parameters are transferable to the lower lying excited states. Modifying the equilibrium bond orders and bond lengths allows the MM calculation to respond to the migration of charge.During the equilibration steps we assume the system to be in its electronic ground state, after which we excite the system to the first excited state and allow the system to respond to the change in the electronic density within the adiabatic Born-Oppenheimer approximation. It is important to note that the excited state we prepare is not the state which carries the most oscillator strength to the ground nor do we account for nonadiabatic surface hopping-type transitions in our approach [31–33]. The dynamics simulations shown reflect the longer-time fate of the lowest lying excited state populations and sample possible configurations that can be accessed by the system. The combination of a classical MD force field with a semiempirical description of a selected few molecules within the system seems to be a suitable compromise between a fully ab initio approach which would be limited to only a few molecules and short simulation times and a fully classical MD description which would neglect any transient changes in the local electronic density [22].


### 4.2. Estimating State to State Rates

We can estimate the state→ state rates using the model by Bittner et al. [23, 24]. Consider a two-state system with coupling λ in which the energy gap Δ(t) fluctuates in time around its average Δ¯. In a two-state basis the Hamiltonian can be written as(15)H=Δt2σz^+λσx^,where σk^ are Pauli matrices. Note that (15) can be transformed such that fluctuations are in the off-diagonal coupling, becoming (16)H=Δ02σz^+δVtσx^,where Δ0=Δ¯+λ and δV¯(t)=0. The fluctuations in the electronic energy levels are attributed to thermal and bond-vibrational motions of polymer chains which can be related to the spectral density, Sω, via(17)V¯2=δV¯2t=∫-∞+∞dω2πSω.Averaging over the environmental noise, we can write the average energy gap as ħΩ¯=Δ02+V¯2 with eigenstates(18)Ψ+=cos⁡θ0+sin⁡θ1Ψ-=-sin⁡θ0+cos⁡θ1,where tan⁡2θ=V¯/Δ0 defines the mixing angle between original kets. Consequently, by analysing energy gap fluctuations, we can obtain an estimate of both the coupling between states as well as transition rates. To estimate the average transition rate between states the equations of motion for the reduced density matrix for a two-level system coupled to a dissipative environment are used.(19)ρ˙11=-iħVρ21-ρ12-1τ1ρ11ρ˙22=iħVρ21-ρ12-1τ2ρ22ρ˙12=-iħVρ22-ρ11-1Tdρ12-Δ0iħρ12ρ˙21=iħVρ22-ρ11-1Tdρ21+Δ0iħρ21.τ1 and τ2 have been introduced as the lifetimes of each state and Td is the decoherence time for the quantum superposition. The decoherence time can be related to the spectral density via Td-1 = V¯/ħ. Taking Td to be short compared to the lifetimes of each state, we can write the population of the initial states as(20)ρ11t=exp⁡-1τ1-k,where k is the average state to state transition rate. If we integrate over all time we obtain an equation of the form(21)∫0∞ρ11tdt=1τ1+k-1suggesting a form for the exact solution of (6). Taking the Laplace transform of the equations of motion (6) and assuming that our initial population is in state 1 (ρ11(0)=1) the equations of motion become a series of algebraic equations(22)-1=-iħVρ21-ρ12-1τ1ρ110=iħVρ21-ρ12-1τ2ρ220=-iħVρ22-ρ11-1Tdρ12-Δ0iħρ120=iħVρ22-ρ11-1Tdρ21+Δ0iħρ21which produces a rate constant of the form(23)k=2V¯2ħ2TdTdΔ0/ħ2+1.


## 4.1. Simulation Methods

Our simulations employ a modified version of the TINKER molecular dynamics (MD) package [25] in which the MM3 [26] intramolecular bonding parameters are allowed to vary with the local π-electronic density as described by a Parisier-Parr-Pople (PPP) semiempirical Hamiltonian [27, 28]. Similar approaches have been described by Lobaugh and Rossky [29] and Muller and Warshel [30] to include electronic dynamics into an otherwise classical force field description.We assume that the nuclear dynamics of all atoms in the system evolve according to a classical description(2)mir¨it=Firt,Pt,where Firt,Pt is the force acting on atom i at time t and Pt is the instantaneous electronic density for the π-electrons. The Fock operator for the π electrons is given by(3)Fμμ=Wμ+12Pμμ∘μμ∣μμ+∑ρ≠μPρρ-ZρVrμρ,Fμν=βμνrμν-12Pμν∘μμ∣νν,where Wμ is the valence state potential for a given atom type and we make the zero-differential overlap (ZDO) assumption to write the electronic repulsion integrals as(4)μν∣ρσ=μμ∣ρρδμνδρσ.The Greek indices correspond to atom-centred basis orbitals. Note that we specifically distinguish P∘ as the ground state (Hartree-Fock) electronic density. Both the hopping integrals βμν and Coulomb repulsion integrals are determined using empirical formulas and depend upon the intermolecular bonding between atoms, accounting for nonplanarity. The molecular orbitals are eigenstates of F(5)Fψn=ϵnψn,where ψn is given as a linear combination of the atom-centred orbitals(6)ψn=∑μcnμχμ.From the cnμ coefficients, we obtain the ground state density matrix(7)Pμν∘=2∑nocccnμcnν.The electronic energy is given by(8)EHF=Tr⁡FP∘and is minimized by imposing the condition(9)F,P∘=0.Electronic π→π∗ excitations are subsequently introduced by using configuration interaction (singles) (CI-S) theory where we write the excited state wave function as a linear combination of elementary excitations between occupied and unoccupied molecular orbitals(10)Ψk=∑h∈occ∑e∈unoccCkehψhψe=∑μν∑ehCkehceμchνχνχμ=∑μνCkμνχνχμ,where the {Ckμν} coefficients give the electron/hole amplitude in the atom-centred basis. These are eigenstates of the CI-S Hamiltonian and we used a total of 10 occupied and 10 unoccupied Hartree-Fock molecular orbitals to construct the electron/hole configurations for the CI calculations. The excited state bond charge density matrix is constructed by assuming that electron densities are added to the unoccupied orbitals and hole densities are subtracted from the Hartree-Fock ground state; namely,(11)Pμν=Pμν∘+Pμνe-Pμνh,where(12)Pμνe=∑σCkμσCkνσPμνh=∑σCkσμCkσν.In the former, we sum over the occupied orbitals (holes) and in the latter we sum over unoccupied orbitals. Pμν is interpreted as the “bond charge” matrix whereby the diagonal elements give the net π-electron count on site μ and the off-diagonal Pμν terms give the π-bond order between sites μ and ν. Consequently, excitations systematically move electron density from atom to atom and change the local bond orders.Within the force field, bond force field parameters are modified by assuming a linear interpolation between single and double bonds. For example, we write the bond force constant between C atomsi and j as(13)k=k2-δk1-Pij,where k2 is the force constant for a C=C double bond when Pij=1, k1 is the bond force constant for a C-C single bond when Pij=0, and δk=(k2-k1). In a similar definition, we interpolate the equilibrium bond lengths between single- and double bonds as(14)l=l2+δl1-Pij.This treatment is not limited to C-C bonds but is applied to all heavy atoms contributing to the π system. While the MM3 parameter set was originally intended to account for the electronic ground state, we assume that the parameters are transferable to the lower lying excited states. Modifying the equilibrium bond orders and bond lengths allows the MM calculation to respond to the migration of charge.During the equilibration steps we assume the system to be in its electronic ground state, after which we excite the system to the first excited state and allow the system to respond to the change in the electronic density within the adiabatic Born-Oppenheimer approximation. It is important to note that the excited state we prepare is not the state which carries the most oscillator strength to the ground nor do we account for nonadiabatic surface hopping-type transitions in our approach [31–33]. The dynamics simulations shown reflect the longer-time fate of the lowest lying excited state populations and sample possible configurations that can be accessed by the system. The combination of a classical MD force field with a semiempirical description of a selected few molecules within the system seems to be a suitable compromise between a fully ab initio approach which would be limited to only a few molecules and short simulation times and a fully classical MD description which would neglect any transient changes in the local electronic density [22].


## 4.2. Estimating State to State Rates

We can estimate the state→ state rates using the model by Bittner et al. [23, 24]. Consider a two-state system with coupling λ in which the energy gap Δ(t) fluctuates in time around its average Δ¯. In a two-state basis the Hamiltonian can be written as(15)H=Δt2σz^+λσx^,where σk^ are Pauli matrices. Note that (15) can be transformed such that fluctuations are in the off-diagonal coupling, becoming (16)H=Δ02σz^+δVtσx^,where Δ0=Δ¯+λ and δV¯(t)=0. The fluctuations in the electronic energy levels are attributed to thermal and bond-vibrational motions of polymer chains which can be related to the spectral density, Sω, via(17)V¯2=δV¯2t=∫-∞+∞dω2πSω.Averaging over the environmental noise, we can write the average energy gap as ħΩ¯=Δ02+V¯2 with eigenstates(18)Ψ+=cos⁡θ0+sin⁡θ1Ψ-=-sin⁡θ0+cos⁡θ1,where tan⁡2θ=V¯/Δ0 defines the mixing angle between original kets. Consequently, by analysing energy gap fluctuations, we can obtain an estimate of both the coupling between states as well as transition rates. To estimate the average transition rate between states the equations of motion for the reduced density matrix for a two-level system coupled to a dissipative environment are used.(19)ρ˙11=-iħVρ21-ρ12-1τ1ρ11ρ˙22=iħVρ21-ρ12-1τ2ρ22ρ˙12=-iħVρ22-ρ11-1Tdρ12-Δ0iħρ12ρ˙21=iħVρ22-ρ11-1Tdρ21+Δ0iħρ21.τ1 and τ2 have been introduced as the lifetimes of each state and Td is the decoherence time for the quantum superposition. The decoherence time can be related to the spectral density via Td-1 = V¯/ħ. Taking Td to be short compared to the lifetimes of each state, we can write the population of the initial states as(20)ρ11t=exp⁡-1τ1-k,where k is the average state to state transition rate. If we integrate over all time we obtain an equation of the form(21)∫0∞ρ11tdt=1τ1+k-1suggesting a form for the exact solution of (6). Taking the Laplace transform of the equations of motion (6) and assuming that our initial population is in state 1 (ρ11(0)=1) the equations of motion become a series of algebraic equations(22)-1=-iħVρ21-ρ12-1τ1ρ110=iħVρ21-ρ12-1τ2ρ220=-iħVρ22-ρ11-1Tdρ12-Δ0iħρ120=iħVρ22-ρ11-1Tdρ21+Δ0iħρ21which produces a rate constant of the form(23)k=2V¯2ħ2TdTdΔ0/ħ2+1.



---

*Source: 6834908-2018-04-10.xml*