Molecular Modeling and Simulation of Conjugated Polymer Oligomers

In accord with earlier studies, the exciton self-traps in the middle of the chain. ... the double and single bonds in the vinylene junctions and p-phe...
0 downloads 0 Views 445KB Size
J. Phys. Chem. B 2008, 112, 4983-4993

4983

Molecular Modeling and Simulation of Conjugated Polymer Oligomers: Ground and Excited State Chain Dynamics of PPV in the Gas Phase Fabio Sterpone* and Peter J. Rossky Department of Chemistry and Biochemistry and Institute for Computational Engineering and Sciences, UniVersity of Texas at Austin, Austin, Texas 78712 ReceiVed: December 18, 2007; In Final Form: February 15, 2008

The ground and excited state dynamics of poly(p-phenylenevinylene) (PPV) chains is studied through an implementation of mixed quantum/classical molecular dynamics simulation. The model used in the simulations combines the semiempirical Pariser-Parr-Pople (PPP) Hamiltonian to treat the π molecular electronic structure with a mechanical force field capturing all other aspects. Nuclear degrees of freedom are treated classically. We first validate the model by simulating PPV chains of various length, and evaluate the absorption spectra. The thermal disorder contribution to the breadth of the first absorption band is estimated to be 0.2 eV at T ) 300 K. To investigate the relationship between the emission and chain conformation, we simulate an isolated ten unit chain of PPV in the ground and the lowest excited state. The emission spectrum, red-shifted with respect to absorption of about 0.2 eV as found in experiments, shows a structured line shape that we relate to the photoinduced CC bond distortions. In accord with earlier studies, the exciton self-traps in the middle of the chain. We introduce two collective variables that reflect geometrical distortion, and find these to be effective in describing the contribution of chain conformation to the emission spectrum. The collective variables are also shown to be effective in describing the bond relaxation dynamics after photoexcitation. Such a relaxation is found to occur in approximately 100 fs and is guided by a compensatory release of energy between the double and single bonds in the vinylene junctions and p-phenyl rings. Finally, we find that the chain has a very slight preference for a more planar conformation in the excited state, compared to the ground state. However, the thermal motions induce the chain to explore out-of-plane conformations in both the ground and the excited states with an amplitude significantly greater than this difference.

I. Introduction Conjugated polymers possess the important optoelectronic properties of semiconductors.1 A new class of devices such as organic light emitting diodes (OLED),2 photovoltaic cells, transistors,3 and flexible displays have been conceived and produced from this class of polymers. However, the general usefulness of such technological devices depends on the competitive cost of engineering and production and on their efficiency. The former aspect is guaranteed by the low cost and flexibility of the plastic material manufacturing, whereas the latter is currently compromised by the consequences of electronic interactions occurring between the polymer chains assembled in the thin films constituting the devices.4 As a result there is great interest in understanding the fundamental interactions and processes occurring in such films. The films of organic conjugated polymers are produced through spin-casting processes from concentrated polymer solution. It has been shown that the optoelectronic properties of the films depend on the precursor solution,5-7 apparently due to aggregation.6 To obtain insight into how the single chain conformation and the relative packing of polymer chains affect the optoelectronic properties, considerable effort has been devoted to the characterization of polymer solutions,8-10 polymer films,5,6 and polymers in glassy matrices11,12 and in confined environments.13-15 As a general feature, it has been demonstrated that interchain interactions such as exciton delocalization * Present Address: CASPUR, Via dei Tizii 6B, 00185, Rome, Italy.

and migration16-18 or polaron pair creation and separation19 play crucial roles in quenching the photoluminescence (PL) properties and lead to a red shift of the emission spectra4 with respect to that of an isolated chain. Furthermore, it is now accepted4 that during the production process, defects occur along the polymer chains (for example cis-conformation, tetrahedral sp3 sites, chain ends, and carbonyl groups) modifying the electronic behavior of polymers in the films. Such defects indeed not only alter the average electronic and mechanical properties of the materials, but also can induce spatial heterogeneity of the photoluminescence activity. To distinctly separate the competitive intra- and interchain interactions, novel experimental approaches have been developed.13,14 Theoretical studies have also contributed greatly to understanding the relationship between optical properties and chain conformation and packing. A first class of studies focused on the casting procedure. For example, Monte Carlo simulations have been used to study the thermodynamics of polymer collapse.20,21 The thermodynamic stability of chain defects induced by casting and their effect on the photophysics has been studied at the atomistic level by the AM1/PPP ab initio method in the case of MEH-PPV and PPV isolated chains.22 It was demonstrated that defects may contribute exciton trapping on shorter segments of the chain. A second class of studies has characterized the excited states of several conjugated oligomer types, monitoring the change in the conformation and the role of packing interactions.23-30 A collective electronic oscillator (CEO) method was used to study the excited state of a single

10.1021/jp711848q CCC: $40.75 © 2008 American Chemical Society Published on Web 04/02/2008

4984 J. Phys. Chem. B, Vol. 112, No. 16, 2008 PPV chain and PPV aggregates.27,28 Most notably, Bre´das and co-workers have carried out a number of studies using different quantum chemical methodologies, providing considerable understanding of chain deformation after photoexcitation and the impact of interchain interactions on charge and energy transfer.23,24 One main focus of these theoretical works is the relationship between spectral properties and chain conformation. Vibrational features of the absorption and emission spectra have been reproduced theoretically by introducing vibronic modes via the Franck Condon progression of vibrational transitions. In the case of PPV oligomers, one finds that at least two effective modes that describe the CC double bond stretching in the vinylene junctions and the CC bond fluctuations in the phenyl rings must be included to reproduce spectra.31,32 A harmonic vibronic analysis of PPV oligomers by Mu¨llen et al.33 showed that low-frequency modes associated with out-of-plane ring libration are primarly responsible for the Stokes shift of the emission spectrum with respect to absorption. Moreover, the asymmetric nature of emission with respect to absorption was traced to mode mixing and frequency shifts between the ground and the excited state. The vibronic fine structure of PPV oligomer emission spectra has been reproduced also by the CEO method.28 Finally, a last class of studies has focused on the natural disorder of the chain conformation. Generally a conjugated polymer chain is idealized as a planar structure. However, steric effects and thermal motion induce out-of-plane conformations. Such distortion from planarity has been monitored, for example, by modeling phenyl ring rotation of conjugated chains in dense polymer films.34 The model was used to describe electrical transport in such films with specific reference to experiments performed on MEH-PPV and PFO films. A random growth algorithm was instead used to simulate the nonradiative excitation energy transfer in PPV chains both in aggregates and linearly confined in nanopores.35 Recently, Molecular Dynamics (MD) simulations with a posteriori analysis of absorption have been performed on isolated PPV chains.36 The classical nuclear motion was simulated by using a mechanical force field and semiempirical electronic structure calculations were carried out on the MD generated configurations. The thermal contribution to the breadth of the absorption has been estimated at various temperatures. It was also shown that the vibronic fine structure underlying the absorption spectra is completely blurred at ambient temperature. A similar a posteriori analysis was used by others to investigate chain properties and oligomer nanoconfinement in the case of PPV and MEH-PPV.37 The effects of thermal disorder on the optoelectronic properties of conjugated polymers, in solution or films, appear both crucial and elusive, and a natural advance in the theoretical investigation of such systems is based on the possibility of combining a realistic description of the electronic structure of the molecules, a model for the molecular flexibility and dynamics, and a reliable description of the environment surrounding such polymers. The mixed quantum/classical molecular dynamics (MQCMD) approach is a powerful tool to address such problems, because it allows one to simulate a quantum subsystem in very complex environments (solution or polymer films) treating the latter mechanically and classically. In this work we present such a model and, as an initial application, the results from MQCMD simulation of isolated PPV chains in the gas phase. The model is based heavily on the seminal work of Warshel and Karplus,38 introduced to treat both the electronic π structure and the nuclear degrees of freedom of oligomer chains. With this model, we are able to

Sterpone and Rossky investigate the role of the chain disorder on the spectroscopic properties as well as the conformational changes induced by photexcitation. In particular, the line shape of the emission spectrum and photoexcitation dynamics is rationalized via two collective variables that naturally reflect the geometrical relaxation induced by photoexcitation. The power of the MQCMD approach is that elaboration of the polymers with pendant groups, the treatment of confining geometries, or solvent effects can be added at little extra computational effort. Extension to nonadiabatic dynamics is also readily accessible.39 In forthcoming works, we will address the issue of interchain interactions and excited state relaxation. This paper is organized as follows: In Section II, we give the details of the model and the simulation protocols. In Section III, we summarize the main results of simulations carried out in the ground state of PPV chains of various length and we present static structural and spectral properties of a ten monomer PPV ((10)PPV) in the lowest excited state. Section IV emphasizes the dynamics of the (10)PPV chain for ground and excited states. Section V presents the conclusions. II. Method The model we assemble to describe the potential energy surface guiding the dynamics of a conjugated polymer chain is discussed in the following. The Hamiltonian is based on the separation of the electronic and nuclear degrees of freedom. In the case of conjugated molecules, the quantum treatment of the electronic structure can be approximately limited to the π electronic system. The σ electrons, lone pair electrons, and inner shell electrons can be combined in effective parameters. Here, the electronic structure is described by the Pariser-Parr-Pople (PPP) Hamiltonian40-42 while the nuclear core interactions with one another are modeled by a mechanical force field. Our Hamiltonian is based on the seminal QCFF/PI method developed by Warshel and co-workers38,43 and the electronic parametrization of the PPP Hamiltonian proposed by Lobaugh and Rossky for the MQCMD simulation of a solvated dye molecule in the ground and excited states.39,44 Our model differs in specific elements from this underlying foundation, and these differences are highlighted in each case below. Here, we will refer to all components of the quantum system environment as “solvent”. For a conjugated molecule in this classical solvent the total potential function VTOT is naturally written as a sum of four contributions:

VTOT ) Vπ + Vc-c + Vs-s + Vs-c

(1)

where Vπ is the electronic part of the potential including the electron-electron and electron-core interactions, as well as interactions of the electrons with the environment when an enviroment is present. The Vc-c term includes all the intramolecular interactions between the nuclear cores such as the electrostatic and the van der Waals interactions as well as the bond, bending and torsional potentials; Vs-s and Vs-c represent the solvent-solvent and solvent-solute interactions, respectively. Precisely because of the quantum nature of the solute, the electrostatic interactions between the solute and the solvent can be divided into (i) the electron-solvent interactions which are included in the Vπ term and (ii) the core-solvent interactions included in the Vs-c term. The Vs-c term includes also the van der Waals interactions between the solvent molecules and the solute atoms. We alert the reader that in the present work the simulations are performed in the gas phase, thus the Vs-s and

Simulation of Conjugated Polymer Oligomers:

J. Phys. Chem. B, Vol. 112, No. 16, 2008 4985

Vs-c potentials are zero and no electron-environment contributions are included in the Vπ term. A. The PPP Hamiltonian. The Vπ term is written as:45

Vπ )

1

∑ Pµν(Hµν + Fµν) 2 µν

(2)

where the summation is over the atomic sites that contribute to the π system. The terms Pµν, Fµν, and Hµν are the bond-order matrix and Fock and one-electron core matrix elements, respectively. The bond-order matrix elements are given by n

Pµν ) 2

∑i ciµciν

(3)

with ciµ and ciν are the molecular orbital coefficients at the sites µ and ν, and the summation runs over all the n occupied molecular orbitals. The PPP method provides the following expression for the Fock matrix elements:46

1 Fµν ) Hµν - Pµνγµν 2 1 Fµµ ) Hµµ + Pµνγµµ + PFFγµF 2 F*µ



(4)

where the elements of the matrix γµν describe the two-electron repulsion. The functional form chosen to describe the two-electron repulsion integral is the Magata-Nishimoto relationship:47

γµν )

e2 Rµν + aµν

aµν )

e2 γµµ + γνν

(5)

where e is the electron charge, Rµν the distance between the atoms µ and ν, and γνν a parameter describing the one-center repulsion. In general this expression has been demonstrated to successfully reproduce molecular spectra.46 The off-diagonal elements, Hµν, of the one-electron core matrix are treated in PPP theory at a comparable level of approximation as Hu¨ckel theory:

Hµν ) βµν

µ, ν covalently bonded

Hµν ) 0

otherwise

(6)

where βµν is the resonance parameter. Thus, only bonded sites contribute to the off-diagonal elements of the H matrix. The parameter βµν can be evaluated for a fixed distance between the π atoms by the Linderberg approximation.48 However, because the bond between two π atoms will be described by a flexible potential,38 the resonance integral is fitted by a numerical expression dependent on the bond distance Rµν (see Appendix E in ref 39):

βµν ) exp(-µµν(Rµν -

µν Req µν))[β1

+

βµν 2 (Rµν

-

Req µν)]

(7)

TABLE 1: Parameters Used in the PPP Hamiltoniana atom



γµµ

βC-X 1

βC-X 2

µC-X

eq RC-X

C

-11.17

9.81

-2.438

-0.9874

2.035

1.397

a The Rµ, γµµ, and are given in eV, eq µC-X is given in Å-1, and RC-X is given in Å.

βC-X 1

βC-X 2

is given in eV/Å,

The term Rµ is a parameter giving the energy of the orbital µ in vacuo in its appropriate valence state. The second term represents the electron core repulsion with ZF, the core charge at the F site. This charge is equal to the number of 2p electrons the atom donates to the π system. In our case, because we consider only singly unsaturated carbon atoms, ZF ) 1. When a molecular environment is present, eq 8 is supplemented by a term accounting for the electron-environment coupling (see eq 9 in ref 44). Within this formulation, the PPP Hamiltonian for conjugated hydrocarbons depends on six types of parameters: the one-center repulsion parameter γµµ, the fitting parameters of the Linderberg µν approximation of the resonance integral (βµν 1 , β2 , µµν, and eq Rµν), and the valence ionization parameter Rµ. The one-center repulsion parameter γµµ as well as the valence ionization parameter Rµ can be deduced from the electronegativity data published by Hinze and Jaffe´;49 γµµ is calculated as the difference between the ionization potential and the electron affinity: γµµ ) I - A.38,50 To coherently match the electronic part with the mechanical force field that describes the intramoµν lecular core-core interactions, the parameters βµν 1 , β2 , µµν, eq and Rµν are derived from the functional results for the same quantities by Warshel and Karplus in ref 38. The set of parameters used in this work refers to the conjugated carbon atoms and are reported in Table 1. B. The Mechanical Force Field. The mechanical interactions are described, except as noted below, by the empirical force field of ref 38. In general, the components of a conjugated system can be separated into saturated atoms, those not contributing to the π electron system, and π atoms that contribute to the electronic structure. In the present PPV oligomers, all the carbon atoms have sp2 hybridization and contribute to the electronic structure. Thus, the Vc-c term includes only the mechanical interactions between π atoms (see eq 34 in ref 38). The potential Vc-c is composed of bonded and nonbonded interactions: B NB Vc-c ) Vc-c + Vc-c

(9)

We first discuss the bonded potential, which includes the potentials for bond stretching, the angle bending, and the torsional rotations: B ) Vc-c

1

M(bi) + ∑ [ki(ai - a0)2 + 2Da] + ∑ 2 i∈CH i∈CC 1

1

∑ Kθ(θi - θ0)2 + 2 i∈torsion ∑angles 2 i∈bond angles 1 2

K(2) φ cos 2φi +

∑ m

1

K(1) φ cos φi +



Kχ(χk - χ0)2 + 2 k Kθθ′(θm - θ0)(θ′ - θ′0) cos φm (10)

The diagonal elements Hµµ are given by the PPP method:

Hµµ ) Rµ -

∑ ZFγµF F*µ

(8)

The first term represents the bond stretching of two π atoms and it has the form of a Morse potential M(bi), which is a physically appropriate description of a bonding potential:

4986 J. Phys. Chem. B, Vol. 112, No. 16, 2008

Sterpone and Rossky

M(b) ) Db[exp{-2R(b - b0)} - 2exp{-R(b - b0)}] (11) with b the bond distance. As the second term in eq 10 shows, the harmonic potential is used for the bonds involving a hydrogen atom. In this case the CH bond distance is referred to as ai. The third term represents the angular bending, with θi being the angles CCC or CCH. The torsions φi around a CC bond are described by the one-fold (K(1)) and two-fold (K(2)) terms. The harmonic term Kχ(χ - χ0)2 accounts for the motion of the out-of-plane angle χ defined by four atoms A, B, and C attached D:

cos χ )

(

)(

)

eAC × eBD eBD × eDC sin θADB sin θBDC

(12)

Finally the last summation considers the two angles θ(XCC’), θ′(CC′Y) of the CC′ bond whose motion is coupled to the torsion φ around the same bond. All the parameters but one entering in the above expressions are taken from the parametrization of ref 38. In ref 38 the core repulsion function γµν as a function of the distance Rµν is given by the Mataga-Nishimoto relationship plus an empirical exponential term not used here. To obtain the correct average bond distance of the CC bond, we modified the equilibrium value b0 of the CC bond; we shift the minimum of the Morse potential of eq 11 from b0 ) 1.46 Å to b0 ) 1.54 Å. This new parametrization allows us to reproduce the average bond distance obtained by the original QCFF/PI force field as is evident in the results (see below). The nonbonded term includes the electrostatic and steric interactions: NB ) Vc-c

γµνZµZν + ∑ AeλR ∑ µν µν

µν

- BRµν-6 + 1-3

1

F(Rµν - R0µν)2 ∑ µν 2

(13)

The electrostatic repulsion between two cores is modeled by a Coulombic model, assigning an integer charge of Zµ at each site µ. However, for consistency with the expression used for the electron-core and electron-electron interactions, the distance dependency is given by the two-electron repulsion matrix elements γµν. Finally, the van der Waals terms include a summation over all the nonbonded distances and a final, harmonic term in which only the 1-3 nonbonded distances are considered. The parameters A, B, λ, F, and R0 51 are the same as in ref 38. To include the effect of the torsional rotation around a given π bond on the bond conjugation, the resonance integral is generally modulated by a factor that varies with the angular rotation. In our model we use the following expression:44

βµν ) β0µν cos θ

(14)

where β0µν is the resonance integral for the planar conformation. For simplicity, the rotation-conjugation coupling is included only for the CC bonds along the vinylene junctions of the PPV chains. C. The Forces. The computation of the forces from our Hamiltonian needed to simulate the dynamics of our systems, either in the ground or in the excited state, poses two problems. First, in the ground state dynamics we need to compute the derivatives of each individual term discussed in the previous paragraphs. All but one, the bond order matrix of eq 3, have a

simple analytical expression that allows efficient computation of the forces. The lack of an analytical expression for the bond order matrix is bypassed by the fact that the molecular orbital (MO) coefficients are calculated through a variational method that minimizes the total electronic energy to obtain the ground state SCF energy, or Vπ.38 Hence, the gradient of the MO coefficients does not contribute to the nuclear forces.38,45 The second task concerns the calculation of the forces in an excited state of the system. In this case the reader is referred to ref 39 where all of the details are provided. Here we only point out that we use configuration interaction with single excitation (SCI) to calculate the excited states. This method is standard for the study of the present systems.23,52,53 Using this method the CI energy is a variational minimum with respect to the CI eigenvector coefficients so that the gradient of the CI coefficients does not contribute to the nuclear forces.54,55 Thus, within this approach, one needs to evaluate only the gradient of the molecular orbital coefficients to obtain the electronic excited state forces acting on the nuclei. The coupled perturbed Hartree-Fock (CPHF) equations56 and the Z-vector method (see Appendix A of ref 39) can be used for this. D. Simulation Protocol. In this initial work, we simulated a set of (n)PPV oligomers of various chain length n in the ground state and a (10)PPV chain in the lowest excited state. The system evolution in all cases is carried out by using the velocity Verlet algorithm with a time step of 1 fs, which provides a relative error in the computation of the total energy of about δE/E = 10-5 in the ground state. In the excited state, because of approximations explained below, the total energy relative error is somewhat larger, δE/E = 10-4. In all the systems, the initial temperature was T ) 300 K and the initial velocities were sampled from a Maxwell distribution. Trajectories are then carried out at constant energy. The trajectories were sampled every 20 fs to compute statistical quantities. For computational reasons we used a finite number of CI states (NCI) to compute the excited state energies and forces. We use NCI ) 50 states in excited state dynamics simulations, and NCI ) 200 in ground state simulations for spectra. It was demonstrated that the absorption spectra were unchanged with an increase in the NCI number. III. Spectra, Structure, and Validation of the Method We validate the method by simulating PPV chains in the ground and excited state and comparison with literature results. The chains were simulated in the gas phase at T ) 300 K. We consider first the ground state properties. PPV chains of various lengths have been simulated in the ground state for 30-40 ps, depending on the chain length. The dependence of the absorption peak on the chain length, the transition dipole orientation, and the CC bond length distribution were investigated. Then, we will turn our attention to the excited state properties. We simulated a chain of (10)PPV in the lowest excited S1 state. To do this, we extracted four independent configurations from the ground state trajectory of the (10)PPV chain and evolved each one of these configurations on the S1 surface for 9 ps. The results we present, photoluminescence spectra and bond length distortion, are obtained averaging over the four independent trajectories. A. Ground State Properties. 1. Saturation Effect. Experimentally, the absorption spectra of PPV oligomers show saturation with increasing oligomer length. Applying the particle in the box model to the conjugated oligomers one expects the absorption energy to decrease linearly as a function of the chain

Simulation of Conjugated Polymer Oligomers:

Figure 1. Absorption spectra for (n)PPV. In the panel labeled (10)PPV, the red and green lines refer to the energy gap e01 computed for two (10)PPV chains, one chain is linear (red) while the second has a cis defect introduced in the middle of the chain. Data are from ref 22; for (3)PPV the red line refers to the computed energy gap e01 of a trans-stilbene molecule in ref 36.

length, as 1/N. It turns out that for lengths larger than a relatively small value, Nc, the spectra peak at approximately the same wavelength. This saturation length Nc is referred to in the literature as the conjugation length of the oligomers. However, it must be pointed out that the microscopic origin of such a critical length is not yet clear. One suggestion is that the limiting conjugation length is an effect of the chain disorder that breaks the planarity of the polymers and then traps the π electronic wave function on a segment of the chain. Others suggest that the conjugation length is a property of the wave function that cannot delocalize over more than a limited number of π bonds, =60,57 because of the electron-phonon scattering.58 While a formal quantum theory of these two views might be phrased differently, it is clear that from a classical nuclei perspective they are not distinct. In Figure 1, we show the absorption spectra computed for PPV chains of various length. We note first that the data show saturation with the increase of the oligomer length. Indeed, for n > 6 the spectra peak approximately at the same energy value (e01 = 390 nm). Thus, we deduce a conjugation length of 6-10 units in agreement with the value of 6-7 units deduced from experiments. Second, our flexible model reproduces very well the energy gap computed by the PPP/SCI Hamiltonian for (10)PPV static structures optimized by the AM1 method (see the green and red vertical lines in the (10)PPV panel).22 For the sake of comparison, we report also in the panel (3)PPV the absorption peak of a trans-stilbene molecule that has been computed by sampling the molecule’s conformation by means of MM3 molecular mechanics and then using (Z)INDO/S-CIS to calculate the vertical electronic transition.36 Our model overestimates the energy gap e01 with respect to experimental results; for each length, e01 is about 0.5-0.6 eV larger than the experimental value. This is probably due to the SCI method, known to overestimate the electronic transitions. In contrast to other theoretical works that use SCI combined with (Z)INDO (e.g., refs 36 and 52), we do not consider here any reparametrization of the method specifically for spectra as was done in the (Z)INDO/S method. This can be done for spectral calculations based on the simulated structures if desired. Figure 2, panel A, shows the offset of our calculations with respect to experiments performed on (n)OPV11 and on (n)t-butyl-PPV59 solutions. Panel B of Figure 2 shows the oscillator strength associated with the S0 f S1 transition as a function of ln(n).

J. Phys. Chem. B, Vol. 112, No. 16, 2008 4987 Data are in very good agreement with other calculations carried out at different theoretical levels.22,36 We note also that the width of the spectrum is reduced with respect to that observed experimentally in dilute solutions for similar conjugated oligomers, but there is no conflict here. The e01 absorption band of our (10)PPV spectrum can be fitted with a Gaussian function of σ = 0.17 eV (a width of 0.3 eV has been used for example for spectral convolution in other works),22 whereas the absorption spectrum of dilute MEH-PPV/chlorobenzene or MEH-PPV/tetrahyrofluoran solutions shows a broad shape of about 0.5 eV width.4 In nondilute solutions and films the spectra are even more broad due to chain/chain interactions and morphology. Though there is not complete agreement,60 the broad absorption spectrum of conjugated polymers seen experimentally is generally rationalized in terms of a superposition of contributions arising from chromophores of different conjugation length (i.e., see refs 61 and 62). In our context, however, this superimposition is absent because the simulated chain is too short to include more than one photoactive fragment. Our spectra are shaped only by the thermal motion of the chains. The contribution of thermal disorder to the spectral breadth has been previously estimated for (n)PPV chains through classical MD simulation to be =0.2 eV,36 very close to our estimate. 2. Transition Dipole. The orientation of the transition dipole is important for understanding the correlation between exciton dynamics and the chain conformation. This orientation relates to the time relaxation of the emission anisotropy following the exciton transfer along the chain.6,63-65 In a PPV chain the vinylene junctions induce an out-of-axis component of the transition dipole µ j so that µ j ‚E h ) cos(θµ), where θµ is different from zero and E h is the longitudinal axis of the chain; see the top of Figure 3. In the literature, the estimate of the angle θµ for the S0 f S1 transition is typically small, θµ ) 5-8°,64-66 with one largest value θµ ) 21°.63 We computed the orientation of the transition dipole with respect to the endto-end distance of the PPV oligomers. We find very good agreement with previous estimates. For static chains (n)PPV, n ) 3, 4, 5, 6, and 10, the angle is θµ = 5°, whereas when the thermal disorder is introduced by the MD simulations, the angle increases slightly up to θµ ) 6.2 ( 2.2°. 3. Bond Length. The distribution of the CC bond lengths is crucial for correctly computing the energy of electronic vertical transitions. The average lengths obtained from our model are d ) 1.405-1.425 Å for the bonds in the rings whereas in the junctions the bond distances are d ) 1.48-1.49 and 1.35-1.36 Å for the single and double bonds, respectively. As discussed in Section II, we used a new parameter for the Morse potential describing the bond stretching between two π carbon atoms. This new parametrization reproduces very well the bond distances computed by simulations by using the original QCFF/ PI force field67 (d ) 1.41-1.42 Å for the bonds in the rings and d ) 1.485 and 1.349 Å for the single and double bonds in the junctions). The bond lengths are also very close to those obtained by classical MD simulations employing the MM3 force field.36 The distribution of the bond lengths is quite uniform along the chain, and only in the external rings deviates from the rest of the chain, a boundary effect. Furthermore, we point out that the bond distribution in the rings slightly deviates from the benzene-like conformation, which is characterized by equivalence of all the bonds of the ring. In our systems, the larger value d = 1.425 Å alternates with the shorter d )1.405 Å (see the scheme in the bottom of Figure 3).

4988 J. Phys. Chem. B, Vol. 112, No. 16, 2008

Sterpone and Rossky

Figure 2. (A) Average energy gap e01 as a function of the inverse of the oligomer length n: black squares from MD simulations, red circles results for (n)OPV oligomers in 2MeTHF solution,11 green circles for (n)t-butyn-PPV in CH2Cl2 solution.59 (B) Oscillator strength as a function of ln(n).

Figure 3. Top: Scheme representing the orientation of the transition dipole with respect to the oligomer axis E. Bottom: The averaged CC bond lengths in the phenyl ring and in the vinylene junction from the ground state simulation. Figure 5. CC bond length distribution along the chain of a (10)PPV monomer from the ground state (black square) and from the first excited state simulations (red circle).

Figure 4. Absorption and photoluminescence spectra computed for a (10)PPV chain. The spectral peak is normalized to 1 to facilitate the comparison.

B. Excited State Properties. 1. Photoluminescence. The photoluminescence spectrum of PPV polymers in solution is red-shifted with respect to the absorption spectrum. Further, it shows a more structured line shape. The emission spectrum is generally traced back to the exciton migration between chromophores:4 the exciton localizing on segments with higher conjugation length. Emission then would arise only from these longer conjugated segments. Conversely, absorption, as discussed in the previous section, is shaped by the contributions arising from segments of various conjugation length. In Figure 4 we compare the photoluminescence and absorption spectra of the (10)PPV oligomer. First, as observed in experiments, the emission peak is red-shifted with respect to the absorption peak. Our spectra are slightly blue-shifted compared to experiments, but the magnitude of the shift between the emission and absorption peaks, =0.2 eV, is comparable to that measured experimentally in dilute solutions of MEH-PPV, 0.2-0.3 eV, although we clearly have no migration.4 The

calculated emission spectrum also shows a noisy but clearly seen shoulder that peaks around 430 nm. A similar shoulder is evident in experimental spectra of several conjugated polymers in solution.4 Such structure is believed to be related to the CC bond stretching. In experiments the shoulder is shifted with respect to the main emission peak by about 0.14-0.15 eV, and in our spectrum the shift is somewhat smaller, closer to 0.1 eV. The position of the shoulder and its relationship with the chain properties will be discussed below. On the basis of these results, we conclude that, although the spectra are blue-shifted with respect to experiments, the molecular structures and the relative position and structure of spectra are in very good agreement with earlier results. 2. Bond Distribution. To understand the photoluminescence properties of PPV chains one needs to characterize the S1 exciton. However, we have not simulated very long chains, and this limitation prevents the study of intrachain migration of the exciton between chromophoric segments. Nevertheless, for (10)PPV, we can observe how the excited state induces changes in the chain conformation and understand exciton localization; the bond length distribution mirrors the exciton localization. In Figure 5 we compare the bond length distribution computed in the excited state to that of the ground state. The distortion of the bond lengths induced by the excitation is localized in the central part of the chain and smoothly damps toward the chain extremes where basically no distortion occurs. This finding, already well-established in previous theoretical work,23,53 highlights the tendency of the exciton to localize in the middle of the chain. The distortion occurs, to different extents, both in the rings and in the vinylene junctions. In the rings, the excited state enhances the difference between the longer (d ) 1.43 Å)

Simulation of Conjugated Polymer Oligomers:

J. Phys. Chem. B, Vol. 112, No. 16, 2008 4989

Figure 6. Bond labels in the phenyl ring and in the vinylene junction.

and the shorter (d ) 1.40 Å) bonds. Thus, the excited state induces a stronger deviation from the benzene-like conformation toward the quinoid-like conformation (for comparison with the ground state, see the scheme in the bottom of Figure 3). On the contrary, along the junction, the excitation tends to reduce the difference between the single bonds at the extremes of the junction (d ) 1.47 Å) and the central double bond (d ) 1.37 Å). The ring and junction distortions can be quantified for each subunit considering the quantities:24,28

1 1 dR ) (dR1 + dR2 + dR3 + dR4 ) - (dR5 + dR6 ) 4 2 1 dJ ) (dJ1 + dJ3) - dJ2 2

observed shoulder at ∼430 nm is related to specific conformational degrees of freedom. For example, the distortion of the junctions can be computed for the subset of configurations having an energy gap e01 smaller than 2.88 eV (i.e., larger than 430 nm). We find that the distortion along the junctions is smaller than average; at the central ring the difference with respect to the ground state increases to ∆dJ ) 0.04 Å. Such a result motivated a more thorough development, discussed next. A. Collective Geometrical Variables. The relationship between the emission spectrum and the chain conformational properties is investigated here by introducing two collective variables that naturally describe the geometrical relaxation occurring after photoexcitation. These collective variables reflect the bond distortion occurring along the chains and involving both the junctions and the rings. We thus introduce the total ring and junction distortions, defined as the sum over all the units of the chain:

(15) 1,NR

(16)

where dRi refers to the bonds in the rings, and dJi refers to the bonds in the junctions (see Figure 6). The distortions dR and dJ are plotted in panels A and B of Figure 7 and compared to the same quantities computed in the ground state. We note that the largest differences between the S1 and S0 states, ∆dJ ) 0.034 Å and ∆dR ) 0.018 Å, occur in the central part of the chain. The magnitudes of such differences are in good agreement with other theoretical calculations;28,53 an AM1 optimized structure of a (4)PPV in the first excited state53 is characterized by a shift with respect to the ground state of ∆dR = 0.015 Å and ∆dJ = 0.04 Å. A similar difference has been computed for the dJ by means of the CEO method.28 IV. Excited State Dynamics The experimental absorption and emission spectra of PPV chains manifest vibrational contributions that mainly originate from CC double bond stretching. Such a vibronic line shape has been reproduced theoretically including vibronic contributions to the spectra, and for PPV oligomers, at least two effective modes (ω ) 1326 and 1667 cm-1) must be included.31,32 Recently, a more complete normal modes analysis was reported by Mu¨llen et al.33 They included up to 20 normal modes and showed that the main vibronic contributions arise from CC double bond stretching modes in junctions and bond fluctuations in rings. Moreover, they showed that mode mixing and a positive frequency shift for the active vibronic modes in the excited state causes the more asymmetric shape of the emission spectrum with respect to absorption. Low-frequency librational modes associated with ring rotation were also indicated as important in excited state relaxation and in inducing the spectral Stokes shift. The vibronic contributions to PL were also investigated by Tretiak et al.28 by their CEO approach and the reported results are in line with the analysis provided by Mu¨llen et al. The vibronic contribution to the absorption of PPV oligomers was also investigated by using (ground state) MD simulations and applying the (Z)INDO method to compute the vertical transitions a posteriori.36 It was reported there that the vibronic fine structure of the absorption spectrum is apparent at low temperature but is lost upon increasing the temperature to 300 K. In the same spirit as these previous works, we have investigated the line shape of the emission spectrum of the (10)PPV chain. In particular, we wanted to clarify whether the

DR )

∑i diR

(17)

1,NJ

DJ )

∑i diJ

(18)

where NR (number of rings) is set to 8, excluding the (undistorted) terminal rings. The number of junctions is NJ ) 9. First, the distributions of the total distortions DR and DJ as well as that of the energy gap e01 are reported in the top and middle parts of Figure 8 for the first excited state and the ground state. In the excited state the distributions are clearly broader than those in the ground state and thus the system explores a larger conformational space in S1. Second, the energy gap surface as a function of DR and DJ was computed and plotted in the bottom part of Figure 8. We plot DR and DJ in the range of 0.15-0.35 and 0.85-1.2 Å, respectively, to focus on the probable regions of DR and DJ distributions. The surface appears as a (somewhat noisy) valley, characterized by a well-defined profile. The isoenergetic surfaces run along the (DR,DJ) diagonal. Hence, one can clearly see that there is a compensatory contribution of the distortions to the energy gap. The isoenergetic curve associated with the shoulder on the long wavelength side of the peak is shown in white, to highlight the contributions from these distortion parameters to the emission structure. B. Distortion Relaxation. The correlation of ring and junction distorsions seen above suggests that they should play a prominent role also in the dynamical relaxation process following photoexcitation. This has been studied here by selecting 10 independent configurations from the (10)PPV ground state dynamics and then simulating the short time relaxation after photoexcitation. The 10 configurations have been simulated for 100 fs in the first excited state and the trajectories have been sampled every 2 fs. The average relaxation of the bond distortions is considered: we computed 〈∆DR(t)〉0 and 〈∆DJ(t)〉0, where ∆Di(t) ) Di(t) - 〈Di〉S1 and i ) R or J, 〈...〉0 means that we average over the 10 initial ground state configurations, and 〈...〉S1 is the average value in the first excited state extracted from the longer time simulations. The averaged trajectory in the (∆DR,∆DJ) plane is plotted in Figure 9. We observe a fast relaxation of the distortions occurring in the first 10 fs and a following clearly oscillatory dynamics in which the system moves along the diagonal around the (0,0) point. The linear correlation between the distortions is associated with a

4990 J. Phys. Chem. B, Vol. 112, No. 16, 2008

Sterpone and Rossky

Figure 7. CC bond distortion computed in the vinylene junctions (panel A) and in the phenyl rings (panel B). Data are presented for the ground state simulation (black squares) and the first excited state simulations (red circles). In panel B we excluded data for the external rings, 1 and 10.

Figure 9. Relaxation of the bond distortion after photoexcitation: averaged trajectory in the distortion space (see eqs 17 and 18). The trajectory is visualized trough a β-spline interpolation. Points are separated by 2 fs. Arrows of different color correspond to trajectory segments of approximately 10-15 fs each.

Figure 8. Probability distribution of DJ, DR (top) and probability distribution of e01 (middle) for the ground state (black line) and the first excited state (red line). In the bottom part of the chart, the energy gap e01 is shown by color scale as a function of the total distortion along the vinylene junctions (DJ) and the phenyl rings DR. See eqs 17 and 18.

Paerson coefficient of r ) 0.78. By inspecting the individual trajectories, we find that at very short time the distortion fluctuations tend to compensate for each other. This corresponds to an alternating transformation of energy from the double to the single bonds in the rings and in the junctions. We point out that within 100 fs, all the bonds basically come to fluctuate around their new equilibrium value.

C. Characteristic Frequencies. The collective coordinates we introduced clearly account for a spread in characteristic frequencies, from vibrational CC stretching to low-frequency collective rearrangements. Here, we discuss the characteristic times of such collective variables. In Figure 10, we present the time correlation functions CJ(t) ) 〈δDJ(t)δDJ(0)〉 and CR ) 〈δDR(t)δDR(0)〉. The correlation functions are computed along the excited state trajectories. We include the short transient segment that follows the photoexcitation, but since the bond relaxation occurs in 100 fs, the correlation function essentially reflects equilibrium sampling. It is evident that CJ(t) initially decays to zero in around 0.6 ps while CR(t) decays in 0.3 ps. Furthermore, in the former case, the decay is underdamped whereas in the latter it is overdamped. The time evolution of the distortion parameters can be further investigated through Fourier Transform (FT) to extract the characteristic frequencies of these degrees of freedom. Both of the distortions show a peak at low frequencies, around ω = 20 cm-1. DJ shows additional peaks between 100 and 400 cm-1 (ω = 100, 310, and 410 cm-1) while DR shows a small peak around 410 cm-1 and an additional structure around 470-480 and 700 cm-1. Hence, these conformational degrees of freedom have distinct characteristic frequencies. Some of these frequencies match quite well those calculated by the AM1 method and measured by inelastic neutron scattering on poly(p-phenylenevinylene) films.68 Comparing our results with those reported

Simulation of Conjugated Polymer Oligomers:

J. Phys. Chem. B, Vol. 112, No. 16, 2008 4991

Figure 10. Panel A: Correlation functions CJ(t) and CR(t) of the bond distortions. Panel B: Fourier transform distortion correlation functions C ˆ J(ω) and C ˆ R(ω).

Figure 11. (A) Velocity-velocity correlation functions. (B) Fourier transform of the velocity-velocity correlation functions. Solid line and dashed line refer to the ground state and lowest excited state, respectively. Inset graph: Fourier transform of the 〈δe01(t)δe01(0)〉.

in ref 68, we can associate the modes of DJ at ω ≈ 100 and 300 cm-1 with the relative displacements of the double bonds in neighboring junctions whereas the modes of DR at 410 cm-1 can be related to out-of-plane motion of the ring atoms and the frequency mode around 700 cm-1 relates to planar vibrations of the ring bonds. Frequencies above 800 cm-1 are not reproduced in these correlation functions because of the sampling frequency, 20 fs. To look at higher frequencies, we computed the velocityvelocity autocorrelation function considering the motion of the carbon atoms. The function is computed from a short supplementary run in the S0 and storing the velocities every 2 fs. The velocity-velocity correlation function is also computed in S1 (data not shown). The function shows, in panel A of Figure 11, a damped decay in the first 0.2 ps, which is characterized by oscillations with a period of T = 0.022 ps. The associated frequency spectrum (panel B) manifests many peaks for ω < 1400 cm-1, a very high and narrow peak around ω = 1570 cm-1 and a secondary peak at 1670 cm-1 that we can associate with the CC double bond stretching and a final peak at higher frequency, ω = 3180 cm-1, which is associated with CH bond stretching. For ω > 3180 cm-1 the spectrum decays continuously to zero. The peak around 1600 cm-1 is related to the oscillatory decay of period T = 0.022 ps. In the S1 state, as a result of the different π conjugation along the chain, the first peak shifts to the value ω = 1600 cm-1 whereas the second, ω = 1670, peak is more pronounced. We recall that in previous theoretical works,53,69 two effective vibronic frequencies have been used to describe the envelope of the emission spectrum of PPV, ω = 1690 and 1290 cm-1, associated to the CC double bond stretch and internal ring mode, respectively. The FT of the energy gap e01 collected from the excited state simulation

(Figure 11B, inset) shows two main peaks at these frequencies. This demonstrates that our model accounts well also for these vibronic contributions, although they are blurred in the spectra by thermal motion. D. Chain Flexibility. The investigation of spectral properties in the case of PPV oligomers deserves a detailed analysis of the out-of-plane motions and specifically of the rotations around the bonds in the junctions. The distribution of the torsional angles around the three bonds of the vinylene junctions have been collected and plotted in Figure 12. The two torsions around the single bonds, which cause the twisting dynamics of two adjacent rings, are denoted by T1 and T2, while the rotation around the double bond of the junctions, which accounts mainly for the bending dynamics of the chain, is denoted by B. In our definition, a planar chain is defined by a value φ ) 180° for each one of these angles. Each distribution has been fitted by a Gaussian function and the fit parameters are reported in the Table 2. After photoexcitation, the twisting rotations, T1 and T2, are associated with a slightly more planar conformation. The average value extracted from the fit actually is φ0 ) 180.7° in the S1 state while it is φ0 ) 181.8° in the S0 state. The T1 and T2 distributions in the excited state are also less broad than those in the ground state, meaning a larger curvature of the torsional potential energy surface as a function of the rotation angle. Those findings correlate qualitatively with reported ab initio calculations of the torsional barrier of a stilbene molecule69 in the ground and first excited state. These calculations have shown that a higher torsional energetic barrier is associated with the excited state. Furthermore, the fluctuations of the rotation around the single bonds agree very well with results from other calculations.36,70 Also in the case of the rotation around the double bonds, we note that the excited state tends to reduce the

4992 J. Phys. Chem. B, Vol. 112, No. 16, 2008

Sterpone and Rossky

Figure 12. Distribution of the torsional angles in the junctions of the (10)PPV. The black line refers to the S0 state and the red to the S1 state. Bold lines indicate the Gaussian fit.

TABLE 2: Gaussian Parameters Obtained Fitting the Distributions of the Torsional Angles in the Junctions of a (10)PPV in the Ground and Excited States system

φ0

σφ

T1 S0 T1 S1 T2 S0 T2 S1 B S0 B S1

181.7 180.7 181.8 180.7 178.1 178.9

11.2 10.9 11.3 11.1 8.8 9.4

out-of-plane deformation of the chain; on average φ0 ) 178.9° in S1 and φ0 ) 178.1° in the S0 state. However, at variance with the rotation around the single bonds, we find that the fluctuations of this rotation are smaller in the ground state than in the excited state. This is an effect of the bond length distortion induced by photoexcitation; the double bond has a larger equilibrium value in S1 than in S0 and the rotation is thus enhanced. The opposite behavior occurs for the rotation around the single bonds, which are shortened after photoexcitation. This analysis confirms qualitatively the results from the CEO method.28 Tretiak and co-workers reported that the PPV chain is more planar in the excited state than in the ground state. A similar picture was proposed in ref 8 to understand the emission spectra of oligo(p-phenyleneethynylene) in solvents of different viscosity. However, it must be pointed out that because of the thermal motion, the chain exhibits considerable flexibility. Actually the end-to-end distance deviates from a planar configuration in both the ground and excited state and relatively large fluctuations (∼10%), mainly caused by the bending dynamics, occur with a characteristic time of about 10 ps. V. Conclusion In this work, we have presented the results from mixed quantum/classical MD simulations performed on PPV oligomers. The simulations are carried out by using a model that consistently describes the π electronic structure and the flexibility of the chain in a framework that is readily extended. We first demonstrated that the model reproduces the main features of the PPV photophysics, including the saturation effect of the absorption spectra, associated with a conjugation length of 6-10 units. Further, the bond lengths and transition dipole orientation agree with previous estimates in the literature.

The excited state dynamics simulation shows the localization of the excitation in the middle of the chain over about 5 units as observed in other theoretical work,23,28 and the bond distortion is in excellent agreement with other theoretical models. For the (10)PPV oligomer, the emission spectrum is red-shifted with respect to the absorption. The emission spectrum computed with only classical dynamics shows structure that is generally related to vibronic contributions. We have introduced a set of collective variables that reflect the geometrical distortion following photoexcitation and permit the characterization of the emission spectrum in terms of the bond length distortions occurring in the vinylene junctions and the rings. We show a clear compensatory contribution of such distortions to the emission spectra. This compensatory contribution also plays a crucial role in the relaxation process following the photoexcitation. We show that an alternating release of energy from the double to the single bonds in the vinylene junctions and rings drives the system during the equilibration in the excited state. We also show that in ∼100 fs all the π bonds have reached a state of oscillation around their average values. The characteristic frequencies of the system motion are computed in the ground and excited state. The bond distortions have characteristic frequencies in the range ω ) 100-700 cm-1. The computed frequencies match those extracted by inelastic neutron scattering experiments, associated with internal modes of the chain (ring and vinylene displacements). Higher frequencies specifically associated with CC stretching are correctly reproduced by our model and contribute to the energy gap fluctuations. Finally, we studied the chain flexibility after photoexcitation. We observed that in the S1 excited state, the rotation around the bonds in the vinylene junctions is slightly more planar than that in the ground state. The fluctuations of such rotations change after photoexcitation as a consequence of the coupling to the bond length distortion. This finding agrees qualitatively with the picture proposed by earlier theoretical and experimental investigations.28,8 However, thermal motion induces the chain to explore very significant out-of-plane configurations in both the ground and the excited states. To conclude, our model permits the exploration of the effects of oligomer chain conformation and motion on the absorption and emission spectra. In particular we have shown that appropriate collective variables that reflect the geometrical distortion in the excited state are effective in describing the emission spectral properties of PPV chains and their relaxation dynamics. This work is a first step toward a more detailed study of conjugated polymers in complex environments such as polymer solutions and films, in addition to the direct simulation of nonadiabatic electronic dynamics. These studies are underway. Acknowledgment. We are very grateful to Prof. A. Warshel and Z. T. Chu for their help in understanding the QCFF/PI and Molaris packages. We thank the NSF for support of this research (CHE-0615173). Partial support by the R. A. Welch Foundation (F-0019) is gratefully aknowledged. References and Notes (1) Friend, R. H.; Gymer, R. W.; Holmes, A. B.; Burroughes, J. H.; Marks, R. N.; Taliani, C.; Bradley, D. D. C.; Santos, D. A. D.; Bre´das, J. L.; Lo¨gdlund, M.; Salaneck, W. R. Nature 1999, 397, 121. (2) Gustafsson, G.; Cao, Y.; Treacy, G. M.; Flavetter, F.; Colinari, N.; Heeger, A. J. Nature 1992, 357, 477. (3) Huitema, H. E. A.; Gelinck, G. H.; van der Putten, P. H.; Kuijik, K. E.; Hart, C. M.; Cantatore, E.; Herwig, P. T.; van Breemen, A. J. J. M.; de Leeuw, D. M. Nature 2001, 414, 599. (4) Schwartz, B. J. Annu. ReV. Phys. Chem. 2003, 54, 141.

Simulation of Conjugated Polymer Oligomers: (5) Nguyen, T.-Q.; Doan, V.; Schwartz, B. J. J. Phys. Chem. B 1999, 110, 4068. (6) Nguyen, T.-Q.; I, M.; Liu, J.; Schwartz, B. J. J. Phys. Chem. B 2000, 104, 237. (7) Huser, T.; Yan, M.; Rothberg, J. Proc. Natl. Acad. Sci. U.S.A. 2000, 97, 11187. (8) Sluch, M. I.; Godt, A.; Bunz, U. H. F.; Berge, M. E. J. Am. Chem. Soc. 2001, 123, 6447-6448. (9) Nesterov, E. E.; Zhu, Z.; Swager, T. M. J. Am. Chem. Soc. 2005, 127, 10083. (10) Grozema, F. C.; Hoofman, R. J. O. M.; Candeias, L. P.; de Haas, M. P.; Warman, J. M.; Siebbeles, L. D. A. J. Phys. Chem. A 2003, 107, 5976. (11) Peeters, E.; Ramos, A. M.; Meskers, S. C.; Janssen, R. A. J. Chem. Phys. 2000, 112, 9445. (12) Hogiu, W.-S.; Peteanu, L. A.; Liu, L. A.; Yaron, D. J.; Wildeman, J. J. Phys. Chem. B 2003, 107, 5133. (13) Cadby, A. J.; Tolbert, S. H. J. Phys. Chem. B 2005, 109, 17879. (14) Nguyen, T.-Q.; Wu, J.; Doan, V.; Schwartz, B. J.; Tolbert, S. H. Science 2000, 288, 652. (15) Bae, S. C.; Lin, Z.; Turner, J.; Granick, S. Polym. Prepr. 2004, 45, 238. (16) Jenekhe, S. A.; Osaheni, J. A. Science 1994, 265, 765. (17) Samuel, I. D. W.; Rumbles, G.; Collison, C. J. Phys. ReV. B 1995, 52, R11753. (18) Jakubiak, R.; Collison, C. J.; Wan, W. C.; Rothberg, L. J.; Hsieh, B. E. J. Phys. Chem. A 1999, 103, 2394. (19) Collison, C. J.; Rothberg, L. J.; Treemaneekarn, V.; Li, Y. Macromolecules 2001, 34, 2346. (20) Doye, J. P. K.; Sear, P. R.; Frenkel, D. J. Chem. Phys. 1998, 108, 2134. (21) Hu, D.; Yu, J.; Wong, K. F.; Bagchi, B.; Rossky, P. J.; Barbara, P. F. Nature 2000, 405, 1030. (22) Wong, K. F.; Skaf, M. S.; Yang, C.-Y.; Rossky, P. J.; Bagchi, B.; Hu, D.; Yu, J.; Barbara, P. F. J. Phys. Chem. B 2001, 105, 6103. (23) Bre´das, J.-L.; Cornil, J.; Beljonne, D.; dos Santos, D. A.; Shuai, Z. Acc. Chem. Res. 1999, 32, 267. (24) Bre´das, J.-L.; Beljonne, D.; Coropceanu, V.; Cornil, J. Chem. ReV. 2004, 104, 4971. (25) Geskin, V. M.; Grozema, F. C.; Siebbeles, L. D. A.; Beljonne, D.; Bre´das, J.-L.; Cornil, J. J. Phys. Chem. B 2005, 109, 20237. (26) Conwell, E. M.; Perlstein, J.; Shaik, S. Phys. ReV. B 1996, 54, R2308. (27) Tretiak, S.; Saxena, A.; Martin, R. L.; Bishop, A. R. J. Phys. Chem. B 2000, 104, 7029. (28) Tretiak, S.; Saxena, A.; Martin, R. L.; Bishop, A. R. Phys. ReV. Lett. 2002, 89, 974021. (29) Ruini, A.; Caldas, M. J.; Busi, G.; Molinari, E. Phys. ReV. Lett. 2002, 88, 2064031. (30) Ferretti, A.; Ruini, A.; Molinari, E. Phys. ReV. Lett. 2003, 90, 864011. (31) Cornil, J.; Beljonne, D.; Shuai, Z.; Hagler, T. W.; Campbell, I. H.; Bradley, D. D. C.; Bre´das, J.-L.; Spangler, C. W.; Mu¨llen, K. Chem. Phys. Lett. 1995, 247, 425. (32) Cornil, J.; Beljonne, D.; Heller, C. M.; Campbell, I. H.; Laurich, B. K.; Smith, D. L.; Bradley, D. D. C.; Mu¨llen, K.; Bre´das, J.-L. Chem. Phys. Lett. 1997, 278, 139. (33) Karabunarliev, S.; Baumgarten, M.; Bittner, E. R.; Mu¨llen, K. J. Chem. Phys. 2000, 111, 11372.

J. Phys. Chem. B, Vol. 112, No. 16, 2008 4993 (34) Yu, Z. G.; Smith, D. L.; Saxena, R. L.; Martin, R. L.; Bishop, A. R. Phys. ReV. Lett. 2000, 84, 721. (35) Claudio, G. C.; Bittner, E. R. J. Phys. Chem. A 2003, 107, 7092. (36) Kwasniewski, S. P.; Franc¸ ois, J. P.; Deleuze, M. S. J. Phys. Chem. A 2003, 107, 5168. (37) Sumpter, B. G.; Kumar, P.; Mehta, A.; Barnes, M. D.; Shelton, W. A.; Harrison, R. J. J. Phys. Chem. B 2005, 109, 7671. (38) Warshel, A.; Karplus, M. J. Am. Chem. Soc. 1972, 94, 5612. (39) Lobaugh, J.; Rossky, P. J. J. Phys. Chem. A 1999, 103, 94329447. (40) Pariser, R.; Parr, R. G. J. Chem. Phys. 1953, 21, 466. (41) Pariser, R.; Parr, R. G. J. Chem. Phys. 1953, 21, 767. (42) Pople, J. A. Trans. Faraday Soc. 1953, 49, 1375. (43) Warshel, A.; Lappicirella, A. J. Am. Chem. Soc. 1981, 103, 46644673. (44) Lobaugh, J.; Rossky, P. J. J. Phys. Chem. A 2000, 104, 899-907. (45) Szabo, A.; Ostlund, N. S. Modern Quantum Chemistry; McGrawHill: New York, 1989. (46) Murrel, J. N.; Harget, A. J. Semiempirical self-consistent molecular orbital theory of molecules; Wiley-Interscience: London, UK, 1972. (47) Nishimoto, K.; Forster, L. S. Theor. Chim. Acta 1966, 4, 155. (48) Linderberg, J. Chem. Phys. Lett. 1967, 1, 39. (49) Hinze, J.; Jaffe´, H. H. J. Am. Chem. Soc. 1962, 84, 540. (50) Bailey, M. L. Theore. Chim. Acta 1969, 13, 56-64. (51) The parameter R0µν is the q0 parameter in eq 34 of ref 38. (52) Cornil, J.; dos Santos, D. A.; Crispin, X.; Silbey, R.; Bre´das, J.-L. J. Am. Chem. Soc. 1998, 120, 1289. (53) Beljonne, D.; Cornil, J.; dos Santos, D. A.; Shuai, Z.; Bre´das, J. L. Excited states in poly(paraphenyleneVinylene) and related oligomers: Theoretical inVestigation of their relation to electrical and optical properties; World Scientific Publ.: Singapore, 1997. (54) Yamaguchi, Y.; Osamura, Y.; Goddard, J. D.; Schaefer, H. F., III A new dimension to quantum chemistry: Analytic deriVatiVe methods in ab initio molecular electronic structure theory; Oxford University Press: Oxford, UK, 1994; p 471. (55) Pulay, P. In Ab initio methods in quantum chemistry-II; Lawley, K. P., Ed.; John Wiley & Sons Ltd: New York, 1987; Vol. 69. (56) Pople, J. A.; Krishnan, R.; Schlegel, H. B.; Binkly, J. S. Int. J. Quantum Chem. Quantum Chem. Symp. 1979, 13, 225. (57) Gubler, U.; Bosshard, C. AdV. Pol. Sci. 2002, 158, 123. (58) Zhang, C.; Braun, D.; Heeger, A. J. J. Appl. Phys. 1993, 73, 5177. (59) Tian, B.; Zerbi, G.; Schenk, R.; Mullen, K. J. Chem. Phys. 1991, 95, 3191. (60) Scholes, G.; Larsen, D. S.; Fleming, G. R. Phys. ReV. B 2000, 61, 13670. (61) Kohler, B. E.; Samuel, I. D. W. J. Chem. Phys. 1995, 103, 6248. (62) Wood, P.; Samuel, I. D. W.; Schrock, R.; Christens, R. L. J. Chem. Phys. 2003, 115, 10955. (63) Hayes, G. R.; Samuel, I. D. W.; Phillips, R. T. Phys. ReV. B 1997, 56, 3838. (64) Hagler, T.; Pakbaz, K.; Heeger, A. J. Phys. ReV. B 1994, 49, 10698. (65) Spano, F. C. J. Chem. Phys. 2001, 114, 5376. (66) Wagersreiter, T.; Mukamel, S. J. Chem. Phys. 1996, 104, 7086. (67) The simulations have been carried out with the Molaris package. (68) Papanek, P.; Fisher, J. E.; Sauvajol, J. L.; Dianoux, A. J.; Mau, G.; Winokur, M. J.; Karasz, F. E. Phys. ReV. B 1994, 50, 15668. (69) Gierschner, J.; Mack, H.-G.; Lu¨er, L.; Olekrug, D. J. Chem. Phys. 2002, 116, 8596. (70) Capaz, R. B.; Caldas, M. J. Phys. ReV. B 2003, 67, 205205-1.