All-Atom Molecular Dynamics Simulation of ... - ACS Publications

Article pubs.acs.org/Macromolecules

All-Atom Molecular Dynamics Simulation of Temperature Effects on the Structural, Thermodynamic, and Packing Properties of the Pure Amorphous and Pure Crystalline Phases of Regioregular P3HT Orestis Alexiadis and Vlasis G. Mavrantzas* Department of Chemical Engineering, University of Patras and FORTH/ICE-HT, Patras, GR 26504, Greece ABSTRACT: Molecular dynamics simulations have been carried out to probe chain self-organization and structure in regioregular poly-3-hexylthiophene (Rr-P3HT). Taking into account recent experimental findings that provide evidence for a semicrystalline ordering in P3HT where crystalline lamellae are periodically separated by interlamellar amorphous zones, we have separately simulated the pure crystalline and the pure amorphous phases of Rr-P3HT, and studied the dependence of their conformational and configurational properties on temperature (from 225 to 600 K). All simulations have been carried out with the very accurate, all-atom Dreiding forcefield. In the pure crystalline phase, the system is found to adopt a noninterdigitated and tilted structure irrespective of temperature and chain molecular weight (Mw) consistent with the results of XRD measurements and other theoretical works. Increasing the temperature did not seem to dramatically affect the interchain distance inside the crystal. Temperature therefore should be expected to have a minor effect on charge hoping in the direction of π−π stacking in P3HT lamellae. As far as the amorphous domains are concerned, our detailed MD simulations have allowed us to study temperature effects on the distribution of dihedral angles along the hexyl chain branches and along the P3HT backbone, as well as on all radial distribution functions. Interestingly enough, as the temperature was lowered below approximately 300 K, the disordered state of P3HT in the MD simulations was found to exhibit signs of crystallization implying a transition from a pure amorphous liquid-like phase to a semicrystalline one. An important outcome of our work is the accumulation of a large number of thoroughly equilibrated atomistic configurations representative of the P3HT pure crystalline and pure amorphous structures, respectively, which can serve as excellent starting points for the execution of larger-scale kinetic Monte Carlo calculations for the computation of charge transport properties based on precomputed values of the hopping rates from site to site in the system to quantify dependence on sample morphology.

I. INTRODUCTION The field of organic electronics is a state-of-art technology, which has drawn considerable attention in designing innovative, convenient, and high-performance electronics.1−3 Breakthrough products employed in this technology include organic lightemitting diodes (OLEDs) used in displays and organic field effect-transistors (OFETs).4,5 Organic photovoltaics and fuel cells also make use of conducting polymers for a number of applications.6−10 One class of materials that are investigated for organic electronic applications is the electronic polymers. Conjugated polymers provide outstanding electrical contact and thermal stability, characteristics that can play an important role for the proper functionality of organic-based electronic circuits.11 These materials exhibit an extended π-conjugated organic backbone that gives rise to some unique opto-electrical semiconducting properties. Polythiophenes (Figure 1) represent one broad family of polymeric semiconductors with such characteristics. These conjugated polymers have a relatively rigid backbone with attached alkyl chains that increase the polymer solubility, and they are generally able to develop © XXXX American Chemical Society

layered crystalline structures with separated main and side chains. Among them, regioregular poly-3-hexylthiophene (RrP3HT) is one of the most widely considered conjugated polymers due to its relatively high carrier mobility (0.1 cm2/V s in the ordered state) coupled with high degree of solubility.12,13 In Rr-P3HT all side chains are pointing in the same direction in a head-to-tale head-to-tale (HT−HT) arrangement. As Figure 1 shows, other arrangements include side chains existing in a head-to-tale head-to-head (HT−HH), tale-to-tale head-to-tale (TT−HT) or tale-to-tale head-to-head (TT−HH) conformation, in which case we use the term regiorandom poly-3hexylthiophene (Rra-P3HT). In P3HT, charge transfer from one chain to another is highly anisotropic as a result of chain packing in the crystalline phase (Figure 2), with charge mobility being faster in the direction of conjugation. Upon chain packing, the molecular ordering in the π−π stacking direction is such that π-orbitals belonging to different chains partially Received: October 23, 2012 Revised: January 26, 2013

A

dx.doi.org/10.1021/ma302211g | Macromolecules XXXX, XXX, XXX−XXX

Macromolecules

Article

Figure 1. Schematic representation of a poly(3-alkylthiophene) (P3AT) unit: (a) regioregular (HT−HT) and (b) regiorandom (TT−HH) arrangements. R represents the side chain, and CR = 6 corresponds to a hexyl branch.

phase transition takes place in the range of temperatures between 225 and 400 K. Recent experimental findings from DSC (differential scanning calorimetry) measurements on P3HT samples of relatively low Mw have revealed the existence of two critical temperatures, one at ∼60 °C related with side chain melting and disordering inside the ordered crystal and another at a higher temperature (between 150 and 250 °C) corresponding to main-chain (thus also to crystal) melting of P3HT.19,24 On a larger scale, Rr-P3HT self-assembles into a wellorganized periodic structure of alternating crystalline domains separated by interlamellar amorphous regions.15,25,26 This periodic microstructure exhibits a hierarchical ordering in three main length scales, namely a π−π stacking distance of 3.8 Å between Rr-P3HT backbones, an interlayer distance of 16 Å between Rr-P3HT backbones separated by alkyl side chains, and a semicrystalline lamellar periodicity of 28 nm.25 Highresolution transmission electron microscopy (HR-TEM) methods have revealed a strong molecular weight dependence of chain packing and semicrystalline structure in oriented films of Rr-P3HT.26 These studies showed that the extension of the crystalline domains tends to decrease with increasing Mw ranging from 20 to 50 nm. Moreover the increase of the total lamellar periodicity for higher molecular weights has been attributed mainly to the increase of the size of the disordered interlamellar regions. The size of the crystalline lamellae, on the other hand, is not substantially affected by Mw. Because of technical limitations of available experimental techniques in determining microscopic structure and charge mobility in semiconducting polymers, it is difficult to elucidate the influence of sample morphology on charge transport. Here, accurate computer simulations and theory can play an important role since they can quantitatively address the relationship between polymer structure and conduction properties thereby also aiding in interpreting the experimental data. Several atomistic and quantum simulations have been conducted with Rr-P3HT to provide insight in the analysis of X-ray diffraction experiments and investigate the structure and charge transfer parameters in crystalline and disordered states.27−39 In detail, quantum mechanical (QM) calculations combined with classical molecular dynamics (MD) simulations have been employed to examine the molecular conformation behavior and properties of charge transport of self-assembled P3HT molecules into a well organized lamellar state. Lan and co-workers30 found that the decrement of regioregularity or/ and the increment of temperature when molecules remain in the ordered lamellae causes a deviation of thiophene rings from planarity which can effect a decrease in the charge carrier mobility along the main-chain direction. However, the direction along the chain backbone (and not along the π−π stacking interchain direction) is still the dominating charge transport route among ordered P3HT regions.

Figure 2. Schematic of the “staggered” configuration proposed for the interlamellar structure of the crystalline phase of Rr-P3HT. The same configuration was adopted for the generation of the initial structures needed in the present MD simulations. Molecules are stacked on top of each other in an orthorhombic unit cell with a = 16.8 Å and b = 7.66 Å but in each “stack” every other molecule is shifted by one thiophene ring along the backbone. The alkyl chains are omitted for clarity.

overlap with each other thereby facilitating interchain charge transfer. In the alkyl chain direction, however, the insulating alkyls inhibit charge transport, thus mobility is the lowest in this direction. The equilibrium crystal structure of regioregular poly-3hexylthiophene and its temperature and molecular weight dependencies have been a matter of debate, with several morphologies proposed over the years based on X-ray and scanning force microscopy (AFM) studies.14−21 Unfortunately, the molecular-scale structure of these systems is not easily resolved by any of these methods due to the lack of long-range morphological ordering stemming mainly from the presence of multiple phases. X-ray measurements in principle are restricted to investigating the crystalline regions of polymer films, while the size of the amorphous zones is essentially large to impair charge transport. The crystalline phase of Rr-P3HT exhibits a lamellar structure where the planar backbone thiophene rings of the molecules are stacked on top of each other in a way that can give rise to two different possible arrangements: The first one, referred to as the “aligned” configuration,22 corresponds to stacked molecules with perfectly aligned thiophene backbones on top of each other whereas the second, the “staggered” one,15 suggests that P3HT molecules are shifted by one thiophene ring along the backbone in the stack (Figure 2). Besides chain stacking, the relative orientation too of molecules belonging to different stacks has not been yet clearly resolved. It has been found that the hexyl side chains and the thiophene rings can be oriented parallel to each other23 or form a zigzag structure.17 Another matter of debate is whether or not a crystal−crystal B


Macromolecules

Article

Calculations based on first-principles density functional theory (DFT) have also been performed to investigate the energetics and stability of different configurations for the RrP3HT crystal28,39 and to determine the influence of backbone tilting29 and torsion angle37 on electronic properties. In particular, the ideal staggered and aligned structures were compared at zero temperature by Maillard and co-workers28 who observed that the most energetically favorable configuration is the staggered structure mainly due to the fact that the alkyl side chains in the staggered arrangement are evenly distributed, thus minimizing the repulsion between them. In addition, Melis and co-workers,39 who found that the π−π interaction is an important contributor to the P3HT selfassembling process, predicted that the aligned crystal structure is only favored in the presence of constraints (e.g., of a substrate); in the absence of constraints, the assembling mechanism leads to a staggered bulk structure in a zigzag configuration of the thiophene rings corresponding to different layers. Recently, there has been an effort to investigate the electronic and atomic structure of the disordered phase of P3HT by combining MD and DFT calculations.35 In view of the existing theoretical and simulation studies on Rr-P3ATs and the nonrealistic (very small) systems addressed by quantum mechanics calculations, the present study aims to shed light on issues that are still on debate as far as the microstructure in semicrystalline Rr-P3HT and its dependence on temperature and polymer chain length are concerned. The present paper is the first step of a more systematic and longterm project aiming at predicting morphology and chain arrangement and structure in poly(alkyl−thiophene) semiconductors at the atomic level. Similar efforts have already appeared in the literature: for example, we can mention the hybrid simulation schemes proposed to investigate realistically large P3HT systems40 and perform large scale studies at the atomic and subatomic level.33−35 Contrary to previous simulation studies on the given polymer that exclusively focus either on the crystalline or the disordered phase of P3HT, in this work we present a thorough comparison of the structural and conformational properties between the two phases for a variety of temperatures. Moreover, issues related to the effect of molecular weight on the morphology of the P3HT crystal are addressed for the first time in this paper. Our long-term objective is to develop a hierarchical modeling approach which will allow us to: (1) check if chains form indeed twodimensional foils or layers or if the much debated presence of interdigitation allows for three-dimensional ordering, as both are suggested by experimental studies; (2) analyze structural properties in the ordered and amorphous regions in these systems (also between the nanodomains, which is known to have a big effect on charge-transfer behavior); (3) understand the effect of regioregularity (stereoregular versus randomregular) on self-assembly and quantify the effect of interchain (lamellae) packing and spacing on charge mobility. We propose a three-step approach as follows: (a) Atomistic MD Study of the Pure Crystalline and Amorphous Phases. By conducting detailed atomistic molecular dynamics (MD) simulations in the isothermal isobaric (NPT) and canonical (NVT) statistical ensemble on relatively large Rr-P3HT systems (∼30000 atoms) with molecular weights corresponding to 20 and 30 thiophene rings (20-3HT and 30-3HT) per chain, our aim in this step is to separately address crystalline and amorphous

phases in these systems. The target behind these simulations is to get fully relaxed structures, representative of the crystalline and amorphous phases of P3HT, with realistic density, atomic packing and chain conformation over a wide range of temperatures. This will give us the first insight into the atomistic structure of the two extreme (pure) phases of P3HT. (b) Atomistic Monte Carlo (MC) Study of the Semicrystalline System. Detailed atomistic simulation of a realistic P3HT sample containing both amorphous and crystalline domains (thus, also the interfacial regions) will be carried out through a Monte Carlo algorithm involving a set of very efficient moves capable of simulating the semicrystalline nature of the system (corresponding to the development of both crystalline and amorphous regions in the material). This step is still underway; it will result in a powerful MC code that will allow us to extend simulations to larger crystal dimensions and longer chain lengths of Rr-P3HT, also to generate model P3HT systems characterized by the presence of P3HT chains in both self-organized structures (lamellae) and amorphous domains. (c) Kinetic MC Simulation of the Charge Transfer Properties in the Crystalline Domain, in the Amorphous Domain, and in the Full Semicrystalline Material. Equilibrated configurations representative of the three structures (pure crystalline, pure amorphous, and semicrystalline) will be used as input in kinetic Monte Carlo simulations to predict their charge carrier mobility at different thermal conditions in an effort to assess the role of morphology on electronic properties. Of course, detailed lower-level quantum mechanical simulations will be needed here to compute the different paths for charge transfer along the chain and from one chain to another, also the corresponding hopping rates. This work is concerned with step (a) of this approach.

II. MOLECULAR MODEL AND COMPUTATIONAL DETAILS The molecular model employed in the present work is based on an implementation of the DREIDING all-atom force-field41 combined with CHARMM potential functions42 to describe the constituent bonded and nonbonded energetic terms: Etotal =

∑

K r(r − req)2 +

bon ds

+

∑

∑

K θ(θ − θeq)2

angles

Kd(1 + cos(nΦ − d))

dihedrals

+

∑ impropers

+

K imp(x − xeq)2 +

⎡⎛ σ ⎞12 ⎛ σ ⎞6 ⎤ 4ε⎢⎜ ⎟ − ⎜ ⎟ ⎥ ⎝r⎠ ⎦ ⎣⎝ r ⎠ nonbonded

∑

Cqiqj εr

(1)

To avoid any confusion, we have to mention that both forcefields (DREIDING and CHARMM) share similar potential functions with small differences (the usage of CHARMM potential function is mainly guided by the fact that in the LAMMPS code the choice of CHARMM notation is widely used for different potential terms). Their main difference is in the treatment of the explicit hydrogen bond term appearing in the DREIDING force-field where it is used to describe the C


Macromolecules

Article

Scienomics.44 Although an empirical method, Gasteiger is very fast and produces accurate results in accordance with other ab initio methods if we take into account the discrepancy among these methods too. For example, we refer to the work of Marcon et al.45,46 who have tested different force-field models for crystalline oligothiophenes. For the case of tetrathiophene (T4), the authors report a comparison of the charges of T4 obtained by an electrostatic potential derived charge (ESP) and a distributed multipole (DMA) analysis. They found that the sulfur atoms are positive with respect to DMA and negative according to ESP. To be on the safe side we tested the set of ESP charges combined with the DREIDING force-field by conducting an MD simulation for the 20-3HT system in the pure crystalline phase at ambient conditions of pressure and temperature. The results proved that the different set of ESP charges had no effect on the structural properties of the crystal. We further note that 1−4 pairwise VDW and Coulomb interactions in the estimation of the system’s potential energy were taken into account in our work with a factor of 1. All classical MD simulations were performed with the LAMMPS47 (version 21 Jul 2010) code. The equations of motion were integrated using the velocity Verlet integrator with a time step whose value was varied between 0.5 and 1.5 fs, depending on the simulation temperature and P3HT phase. Long-range electrostatic interactions were calculated by implementing the particle−particle, particle−mesh (P3 M)48 method. The 12−6 Lennard-Jones (describing the VDW interactions) and the Coulombic terms were computed with an additional switching function Sw(r) that ramps the energy and force smoothly to zero between an inner (8 Å) and an outer (13 Å) cutoff:

interactions of a hydrogen atom with electronegative atoms (like N, O, F); of course, this is not the case here (and therefore was not taken into account). In detail, the form of the potential energy function employed in our MD simulations includes four types of bonded-energy contributions (bond stretching, bond angle bending, torsional, and improper dihedrals) and two nonbonded terms (van der Waals (VDW) and Coulomb). The parameters appearing in the bonded terms (Kr, req, Kθ, θeq, Kd, n, d, Kimp, and xeq) and in the van der Waals expression (ε, σ) have been taken from the DREIDING model potential. All the force-field parameters are given analytically in Table 1. System atomic charges (q) were calculated using the Table 1. DREIDING Force-Field Parameters Governing CHARMM Potential Functions of Eq 1 atom types CR CB S H pair coeffsa

angle coeffs

θeq (deg)

carbon in the ring X−CR−X 120.00 carbon in the hexyl CR−S−CR 92.40 branch 109.47 sulfur in the ring X−CB−X hydrogen ε Kd σ (Å) (kcal) dihedral coeffs (kcal)

CR−CR CB−CB S−S H−H

3.473 3.473 3.590 2.846

bond coeffs

req (Å)

Kr (kcal/Å2)

CR−CR CR−CB CR−S

1.33 1.43 1.70

700 350 350

CR−H

0.99

350

CB−CB CB−H

1.53 1.09

350 350

0.095 0.095 0.344 0.015

X−CR−CR−X X−CR−S−X X−CR−CB−X X−CB−CB−X H−CB−CB−H

5.625 0.1667 0.0833 0.1111 0.0000 xeq improper coeffs (deg)

CR−H−CR−CR CR−H−S−CR CR−CR−CR− CR CR−CB−CR− CR

0.0

KΘ (kcal/rad2) 50 50 50

n

d (deg)

2 2 6 3 0

−180.0 −180.0 −180.0 0.0 0.0 Kimp (kcal/rad2) 20.0

Sw(r ) =

(rout 2 − r 2)2 (rout 2 + 2r 2 − rin 2) (rout 2 − rin 2)3

(2)

MD simulations concerning the crystalline phase of RrP3HT were performed in the isothermal isobaric (NPT) statistical ensemble. The Nosé−Hoover thermostat/barostat with corresponding relaxation times equal to 10 and 350 fs were used to control the simulation temperature and pressure. The simulation cell was subjected to periodic boundary conditions along all three directions (XYZ). For the corresponding simulations concerning the amorphous phase of Rr-P3HT, and in order to get fully relaxed and uncorrelated configurations, the crystal was initially annealed at a high

a Cross terms were calculated using the Lorentz−Berthelot (arithmetic) mixing rules.

Gasteiger method43 of partial equalization of orbital electronegativity through the MAPS simulation platform of

Figure 3. “Staggered” initial configuration of the 20-3HT system used in the MD simulations of the crystalline phase. Chain backbones are directed along the X direction. The system consists of 16 sheets distributed in the Y direction (for clarity only four sheets are depicted); each sheet contains four chains of 20-3HT distributed along the Z direction (for clarity, only two chains per sheet are depicted). Values of the lattice parameters are reported in Figure 2. Hydrogen atoms are omitted for clarity. D


Macromolecules

Article

Figure 4. Evolution of the P3HT crystal’s (a) potential energy and (b) density in the course of the MD simulation at 300 and 370 K. (c) Density of the 20-3HT crystal averaged over several MD equilibrated frames as a function of simulation temperature. (d) Evolution of the instantaneous value of the mean square end-to-end distance of the 20-3HT polymer backbone in the course of the MD run, as a function of simulation temperature.

cell as shown in Figure 2. By considering 4 × 16 chains in the Z and Y directions respectively for the 20-3HT system (and 4 × 12 for the 30-3HT system) an initial configuration of 64 chains for the 20-3HT (and of 48 chains for the 30-3HT) ring system was obtained as depicted in Figure 3.

enough temperature (close to 1500 K, where complete equilibration was considerably easier to achieve within the available CPU time) in the isothermal−isobaric (NPT) ensemble and then successive coolings of carefully selected, totally uncorrelated melt configurations in the NPT ensemble took place down to the temperatures of interest. Two different P3HT systems were simulated: the first consisted of 64 chains each one made up of 20 thiophene rings (we call this the 20-3HT system) and the second of 48 chains each one containing 30 thiophene rings (we call it the 30-3HT system). The selection of the number of thiophene rings (20 and 30) along the chain was made having in mind earlier experimental findings,26 according to which the size of the crystalline regions in P3HT along the chain axis is around 8 nm independent of the Mw of the system; from this chain length and above, P3HT chains start to twist and fold, forming amorphous interlamellar zones. On the basis of these reports, in our case we chose one system that is purely crystalline (the maximum length of a 20-3HT chain is below 8 nm) and one that is expected to display chain twisting after a certain number of thiophene rings along the backbone is exceeded (the maximum length of a 30-3HT chain is approximately 11 nm). For both systems, the initial configuration was constructed with all alkyl−thiophene chains in the staggered15 arrangement with noninterdigitated side chains, aligned in an orthorhombic unit

III. RESULTS Crystalline Phase. Exhaustive MD simulations were conducted in the NPT ensemble for the 20-3HT crystal above and below the experimentally observed and reported crystal−crystal transition temperature of ∼60 °C19 (four temperatures were studied: 225, 300, 335, and 370 K) for more than 25 ns using a time step equal to 1.5 fs. Because the initial setup of the crystal was not energetically favorable and far from its equilibrium state (and in order to deal with possible trapping in states of deep local energy minima especially in the lowest simulation temperatures), the initial configuration was first annealed to 370 K until a fully equilibrated structure was obtained (in the sense that the instantaneous values of all structural and energetic quantities were observed to fluctuate constantly around a steady mean value, as it is depicted in parts a, b, and d of Figure 4 for the system’s potential energy, density, and mean square end-to-end backbone distance, respectively). The computed density of the crystalline structure with respect to simulation temperature which is displayed in Figure 4c is E


Macromolecules

Article

directly comparable to experimental findings14,15,49 which report a value of 1.1 g/cm3 at 20 °C for high Mw samples (an order of magnitude higher than the simulated systems). However, the P3HT samples employed in the experimental studies are not perfect crystals but contain amorphous zones too. Thus, the measured experimental density is expected to be smaller than the one corresponding to the pure crystalline phase. This points out to a very favorable comparison with the predictions of the present simulations. The relaxed configuration of the crystal at the end of the MD simulation at 370 K was then cooled down to the lower temperatures of interest. This resulted in fully relaxed structures for the 20-3HT crystal with realistic atomic packing, chain conformational features and potential energy values at all temperatures (Figure 5). Irrespective of temperature, P3HT

of crystal structure on the lattice parameter a (see Figure 2). To this, we conducted short simulations (on the order of a few hundreds of picoseconds) at a low temperature (150 K) starting from initial configurations with three different values of a (two smaller and one larger than the 16.8 Å used here, namely 13, 16, and 18 Å). Interestingly enough, we observe in Figure 7 that even at such a low temperature the crystal structure is destroyed when side chains are initially forced to be strongly interdigitated (a = 13 Å, Figure 7b); even the π−π registry between rings which is the main contributor to the P3HT self-assembly is lost. The same result, although less pronounced, appears for a = 16 Å (Figure 7d). For larger values of a (e.g., 18 Å, Figure 7f), the crystal retains its organized structure resembling the noninterdigitated final configurations of Figure 5. Small deviations in the conformation of side chains in Figure 7f from the “zigzag” structure of Figure 5b should be attributed to the low simulation temperature (150 K) which constrains hexyl branches in the all-trans conformation. The well-ordered structure and organization of the 20-3HT molecules inside the crystalline phase was quantified by analyzing the distribution of a number of characteristic angles and the profile of several pair distribution functions; wherever possible, we also compared against available experimental data and other simulation results. For instance, to interpret the “zigzag” ordering inside the crystal we investigated the distribution of the tilt angle (see Figure 6a) representative of the inclination of the ring planes along the chain with respect to the Z axis; the results are shown in Figure 6b. It is obvious from the plots that there is a preferential value of the tilt angle around the 20° value, which is consistent with the findings of recent computations based on density functional theory.29,30 In Figure 8 we present the dihedral angle distributions for the torsional angles along the side hexyl-chains and in the thiophene rings. In parts b and c of Figure 8, in particular, we report the distribution of the dihedrals corresponding to the first two bonds of the hexyl branches. Irrespective of the temperature, the first bond is observed to be predominantly in the all-trans conformation reflecting the coplanarity with the corresponding 5-mer ring. The rotation of the second bond with respect to the ring plane (by ∼90°) reflects the “zigzag” structure also seen in Figure 5b and Figure 6a. The distribution

Figure 5. Final configuration of the 20-3HT system acquired at the end of a 27 ns-long MD simulation at 300 K: (a) projection on the XZ plane and (b) projection on the YZ plane of the coordinate system. For clarity, hydrogen atoms are not included in the snapshots and only a part of chain backbones is displayed.

chains were found to display a “zigzag”39 noninterdigitated and tilted smectic-like structure between thiophene rings and alkyl branches as shown in Figure 5(b) and clearly demonstrated in Figure 6a. Our result concerning the noninterdigitation of side chains agrees with early but also with more recent diffraction studies of P3ATs17,50−53 as well as theoretical arguments30 and simulation results.32,39 The noninterdigitation of the alkyl branches was further supported by looking at the dependence

Figure 6. (a) Representation of the “zigzag” structure adopted by the 20-3HT chains inside the crystal at the end of the MD simulation at 300 K. The chains depicted have been selected from the YZ projection of Figure 5b; the utilt vector defines the tilt angle formed between the Z axis and the normal to the ring plane. (b) Distribution of the tilt angle for the 20-3HT ring planes at three different simulation temperatures. F


Macromolecules

Article

temperature range. There is, however, a different path to study possible intrachain crystalline phase transition, by looking at the dynamics of the side chains. To this, we have studied the decay of the time autocorrelation function ⟨uBr(t)·uBr(0)⟩ of the unit vector uBr directed from one end to the other end of the hexyl branch and how this accelerates with increasing temperature. Our simulation predictions are shown in Figure 10 from which the following observations can be made: there is practically no relaxation of the side chains for temperatures equal or below T = 335K as the corresponding curves do not seem to drop to zero at all. In contrast, at T = 370 K, this function does go to zero albeit quite slowly. The corresponding relaxation time can be estimated by fitting the simulation curve with a stretched exponential (or KWW) function of the form β ⎛ t ⎞ ⟨uBr(t ) ·uBr(0)⟩ = A exp⎜ − ⎟ ⎝ τKWW ⎠

(3)

and by calculating the integral below the curve. In eq 3, τKWW and β are the characteristic relaxation time and stretching exponent parameters, respectively, and A the amplitude introduced to account for the fast decay of ⟨uBr(t)·uBr(0)⟩ at subpicosecond time scales. The integral below the ⟨uBr(t)·uBr(0)⟩-vs-t curve defines the total correlation time τc characterizing the relaxation of ⟨uBr(t)·uBr(0)⟩ and is given by the following expression:

( )τ

Γ τc = A

1 β

β

KWW

(4)

The value of the relaxation time thus obtained for ⟨uBr(t)·uBr(0)⟩ at T = 370 K comes out to be ∼12 s. Our MD data indicate, therefore, a kind of transition at a temperature between 335 K and 370 K from a state of practically fully frozen (crystalline) side groups to a state of mobile (molten) side groups, which is consistent with the DSC measurements. Additional quantitative information about the microstructure (atomic packing and thiophene ring arrangement) in the P3HT crystal has been acquired from the investigation of the radial distribution functions (RDFs) of selected atom pairs. Figure 11 shows the partial RDFs for the intermolecular (i.e., belonging to different chains) part of the ring center-of-mass, of the sulfur−sulfur, and of the carbon−carbon pairs, and of the intrachain part of sulfur−sulfur pairs, at the two lowest simulation temperatures (225 and 300 K). Both interchain S−S and C−C RDFs exhibit strong peaks (indicative of the well-ordered lamellar nanostructure of P3HT), which as expected become less intense at higher distances or as the temperature is increased. These peaks are more pronounced in the RDF of the ring center-of-mass shown in Figure 11a and provide a better picture of the interchain spacing due to π−π stacking (the b lattice parameter in Figure 2). In Figure 11a, this corresponds to the third peak (the one around 7.62 Å) and reflects the distance from the next nearest center-of-mass neighbor; alternatively, this can be estimated from the nearest neighboring center-of-mass as twice the value of the first peak which is 3.8 Å. The intrachain spacing along the P3HT backbone is quantified by the RDF of sulfur−sulfur pairs along the same chain (intramolecular part), as depicted in Figure 11d. The characteristic distance between every other sulfur atom along the backbone of a P3HT chain corresponds to the second peak and suggests a value of 7.65 Å.

Figure 7. Initial and final (acquired at the end of short MD simulations at 150 K) configurations of the 20-3HT system for different values of the lattice parameter a: (a and b) a = 13 Å; (c and d) a = 16 Å; (e and f) a = 18 Å.

profile of the dihedral angle between successive thiophene rings along the chain at all simulated temperatures is depicted in Figure 8d. Our results demonstrate that, with increasing temperature, the distribution of the rotational angle between the backbone thiophene rings tends to broaden up. In addition, we observe that irrespective of the simulation temperature a significant percentage of the inter-ring dihedral population acquires values which deviate from coplanarity (180°) by about ±20°.30 In the literature, experimental data from DSC measurements on P3HT samples suggest an intracrystalline phase transition due to side chain melting (which causes significant disordering of the P3HT crystal) at around 60 °C. In our MD simulations it is obvious that the structure and dynamics of the branches are governed by the rotation of the second bond which adopts always the gauche conformation irrespective of temperature, as can be clearly seen in Figure 8. The rest of the bonds along the hexyl branch seem to display an oscillation in the fraction of gauche defects between a maximum and a minimum value for every other bond as it is analyzed quantitatively in Figure 9. According to our MD simulation data, therefore, alkyl chains are not predominantly in the all-trans conformation but exhibit a significant percentage of gauche defects even at temperatures as low as T = 225 K. On the basis of the dependency of branch gauche defects on temperature, our simulations do not reveal any signs of imminent status change that would allow us to deduce sound conclusions for a phase transition in the expected G


Macromolecules

Article

Figure 8. Dihedral angle distributions along the hexyl-chain branches and the backbone for the temperatures of 225, 300, and 370 K: (a) description of the dihedral angles, (b, c) first (α-CH2) and second (β-CH2) bond dihedral distributions of the hexyl branches, and (d) inter-ring dihedral (SCCS) distribution among consecutive rings along the P3HT backbone.

Figure 9. Percentage of bonds in the gauche configuration for dihedral angles along the hexyl branch, and its temperature dependence.

Figure 10. Decay of the time autocorrelation function ⟨uBr(t)·uBr(0)⟩ for the unit vector uBr directed from one end of the hexyl branch to the other for the 20-3HT crystal, as a function of simulation temperature.

An investigation of the effect of temperature on the microstructure is attempted in Figure 11e, where the RDF of the ring center-of-mass is presented for the lowest and highest simulation temperatures addressed here. With increasing temperature, there is an obvious shift of all characteristic peaks toward larger distances followed by a drop in their intensity; nonetheless, the increase in the interchain spacing

with temperature is negligible to the degree that it prevents us from arguing that it can induce a decrease in the charge mobility between neighboring chains (for example, the first peak which corresponds to the π−π stacking spacing shifts only H


Macromolecules

Article

Figure 11. Radial distribution functions at T = 225 and 300 K: (a) between the centers-of-mass of rings belonging to different chains (intermolecular part), (b) between intermolecular sulfur pairs, (c) between intermolecular carbon pairs, and (d) between intramolecular sulfur pairs along the chain’s backbone. In part e, the same curves are shown as in part a but for T = 225K and 370 K.

where the final snapshots from the simulations with both crystals (at 370 K) are depicted. In Figure 13a, the effect of molecular weight on the crystal’s instantaneous density is presented; interestingly enough, we observe that the 20-3HT crystal is denser than the 30-3HT one (0.966 g/cm3 vs 0.946 g/ cm3 on average) reflecting a higher degree of packing in the case of the shorter system. This can be explained by looking at the XY projections of the final configurations of the 20-3HT and 30-3HT crystals, which are shown in Figure 12, parts c and

from 3.8 to3.9 Å when the temperature is increased from 225 to 370 K). To obtain additional insight into the microstructure of the P3HT crystal, we conducted MD simulations of a system with a relatively higher molecular weight at T = 370 K. It is the 303HT system consisting of 48 chains each one of which is made up of 30 thiophene rings. This crystal is less extended in the π−π stacking direction (Y direction in the simulation box) compared to the 20-3HT system as can be seen in Figure 12 I


Macromolecules

Article

d. We observe a chain twisting in the case of the longer-chain crystal at distances along the chain greater approximately than 8.0 nm from its one end, a result consistent with experimentally reported observations.26 As a consequence of chain twisting, the 30-3HT system occupies proportionally a larger space volume than the 20-3HT one thereby appearing less dense. We tried to build up on this significant result by attempting to simulate an even longer sample, namely a 48-chain 40-3HT system (P3HT chains made up of 40 thiophene rings) at a very low temperature (150 K). Clearly, the MD method becomes inefficient (in terms of computational time) when large systems with ordered structures are to be simulated, especially at very low temperatures as the one considered here; nonetheless, we got a more pronounced chain twisting for the 40-3HT system as it is clearly depicted in Figure 14; the importance of this result lies in the fact that chain twisting is observed to occur again at the same distance from the chain start (≈8.2 nm) as for the 30-3HT one (≈8.0 nm). The difference in chain packing between the 20-3HT and the 30-3HT crystals is also evident in the final snapshots reported in Figure 12b, where a degree of deviation of the backbone ring rotation with respect to the Z axis (tilt angle) away from 20° is observed; this defect is more pronounced for the chains belonging to the first and last of the four 12-chain layers. We have quantified this result in Figure 13b, where we show the distributions of the tilt angle corresponding to the two different chain length systems: a systematic tilting of the 30-3HT ring planes at angles less than 20° is observed accompanied by a relatively smaller density. Amorphous Phase. For the studies of the amorphous domains, we started the MD simulations from a fully relaxed and uncorrelated configuration with no indications of crystalline remnants, obtained by annealing the equilibrated configuration of the crystal acquired at the end of the NPT MD simulation at 370 K to a very high temperature (1500 K) to ensure perfect crystal melting. This is verified by the graphs shown in Figure 15 displaying (a) the relaxation of the time autocorrelation function for the chain end-to-end unit vector ⟨u(t)·u(0)⟩, (b) the evolution of the instantaneous value of the mean-square chain end-to-end distance ⟨R2⟩, and (c) the corresponding time evolution of the mass density ρ of the

Figure 12. Typical configurations of (a) the 20-3HT system (64 chains) and (b) the 30-3HT system (48 chains) at the end of the MD simulation at 370 K (YZ projections). (c) and (d) XY projections for the 20-3HT and 30-3HT crystals (for clarity only one out of the four layer is depicted). In part d, the longer molecular length causes “breaking” of the crystal and subsequent chain twisting. Hydrogen atoms have been omitted from the snapshots.

Figure 13. Comparing the instantaneous density values in the course of the MD simulation (a) and the average tilt angle distributions (b), between the 20-3HT and the 30-3HT crystals at T = 370 K. J


Macromolecules

Article

Figure 14. Typical configuration of the 40-3HT system acquired after 6 ns of an MD simulation at 150 K: (a) Projection on the XY plane. “Breaking” of the crystal and subsequent chain twisting occurs after a certain molecular length is exceeded. (b) Projection on the YZ plane of the coordinate system. For clarity we have omitted some chains in both layers. Hydrogen atoms are not included in the snapshots.

system. ⟨u(t)·u(0)⟩ is seen to quickly drop to zero within 3.5 ns, demonstrating complete relaxation of the conformational features of the chains associated with their longest length scale. ⟨R2⟩, on the other hand, is seen to exhibit significant fluctuations (±313 Å2) in the course of the MD run, indicative of ergodic and robust sampling by the MD method. As far as the density is concerned, this is seen to reach almost instantaneously a relatively low and steady-state value characteristic of the molten system at 1500 K. After the annealing process is finished, a configuration devoid of any signs of crystal ordering is produced as it can be seen in Figure 16a. A set of different configurations from different trajectory points in the equilibrated regime (see Figure 15, parts a and b) like the one depicted in Figure 16a served as starting points for executing several gradual coolings in the NPT ensemble followed by alternating NVT relaxations in a range of temperatures from 1000 to 600 K and at a pressure P equal to 1 atm, typically for 3−6 ns. Further equilibration at the lower temperatures of interest (600, 500, and 300 K) was in order to get a fully relaxed final structure of the amorphous phase characterized by realistic density values and atomic packing features (Table 2 displays technical information for each simulation temperature concerning the total MD execution time and the corresponding time step). Such a structure is representative of the 20-3HT system configuration at 300 K and is shown in Figure 16b. It is obvious that at 300 K the initial amorphous structure of the 20-3HT polymerwhich was obtained from the MD run at 1500 K; see Figure 16a starts to exhibit signs of crystallization resembling more its semicrystalline nature. This transition is predominantly

reflected in the sudden drop of density; depicted in Figure 17 is the equilibration of the 20-3HT amorphous system with respect to the system’s instantaneous density (Figure 17a) during the gradual cooling to T = 800 K, then to T = 600 K, then to T = 500 K, and finally to T = 300 K. Focusing especially in the curve shown in Figure 17b, we observe an almost linear increase of the system’s density as the simulation temperature is decreased from 800 to 300 K. As expected, the density of the semicrystalline 20-3HT system is significantly smaller than that of the 20-3HT pure crystal at 300 K (0.9675 ± 0.007 g/cm3 versus 1.004 ± 0.01 g/cm3). The predicted density of the amorphous phase at 300 K is in good agreement with the measured monomer 3HT density (0.936 g/mL)54 at the same thermodynamic conditions. Although the available experimental data14,49 report a density for the semicrystalline P3HT system around 1.1 g/cm3 (i.e., a value larger than the one predicted from our pure amorphous phase simulations), we must not forget that a realistic P3HT sample has a much higher degree of crystallinity compared to the simulated high disordered P3HT system. Figure 18 shows the dihedral angle distributions for selected bonds along the hexyl branches and for the inter-ring torsional angle along the P3HT backbone as a function of simulation temperature. We find that overall the hexyl branches retain their structural and conformational features either being in the pure amorphous or in the pure crystalline phase of the P3HT system. For example, from Figure 18a, we observe that the first bond exhibits the same degree of coplanarity with the corresponding 5-mer ring as in the 20-3HT crystal and that the rotation of the second bond with respect to the ring’s plane K


Macromolecules

Article

Figure 15. (a) Decay of the time autocorrelation function ⟨u(t)·u(0)⟩ for the unit vector u directed from one end of the P3HT to the other for the 20-3HT melt, as obtained from the NPT MD simulations at T = 1500 K and P = 1 atm. (b) Time evolution of the mean square chain end-to-end distance ⟨R2⟩ of the 20-3HT model system during the 6 ns-long NPT MD simulation at T = 1500 K. (c) Time evolution of the corresponding density ρ of the 20-3HT melt model system during the 6 ns-long NPT MD simulation at T = 1500 K. Arrows in parts a and b denote frames from the trajectory taken for cooling to lower temperatures.

Table 2. Technical Details for All Simulations Conducted for the Pure Amorphous Phase of P3HT temperature (K) (crystal heating) 1500 1500−1300a 1300 1300−1000a 1000 1000−800a 800 800−600a 600 600−500a 500 500−300a 300

simulation time (ns)

statistical ensemble

time step (fs)

6

NPT

0.5

3 1 4.5 1 2 16 3 11 2 13 3 15

NPT NVT NPT NVT NPT NPT NPT NPT NPT NPT NPT NPT

0.5 0.5 0.5 1 1

Figure 16. (a) Typical configuration of the 20-3HT amorphous system after annealing the crystal at a high temperature (1500 K) for about 6 ns. (b) Semicrystalline structure of the 20-3HT polymer at 300 K acquired at the end of the gradual cooling process that was started from the configuration of part a. For clarity the hydrogens are not included in the snapshots.

a

(Figure 18b) keeps fluctuating around the value of 90° exactly as it does in the crystal case (the only difference between the two phases as far as the conformation of dihedral angles along branches is concerned being that the one corresponding to the

amorphous case is somewhat broader). As expected, the transstate population of the dihedral angle corresponding to the free end of the hexyl branch (Figure 18c) is decreasing with increasing temperature reflecting a more flexible bond at the L

1 1

Gradual coolings.


Macromolecules

Article

Figure 17. (a) Evolution of the amorphous system’s density in the course of the MD simulation at the temperatures of 800, 600, 500, and 300 K, starting from the fully relaxed and low-density configuration obtained at the end of the thermal annealing process at 1500 K (see Figure 15c). (b) Density of the 20-3HT system (averaged over a large number of MD equilibrated frames) as a function of simulation temperature.

Figure 18. Dihedral angle distributions of the amorphous 20-3HT system along the hexyl-chain branches and the ring’s backbone for the temperatures of 300, 500, 600, and 800 K; (a−c) first (α-CH2), second (β-CH2), and last (CH3) bond dihedral distributions of the hexyl branches and (d) inter-ring dihedral (S−C−C−S) distribution among consecutive rings in the backbone.

(±360°) for the inter-ring dihedrals. This latter feature

highest temperatures. In Figure 18d, we display the distribution profile of the dihedral angle along the backbone of the thiophene rings for all simulated temperatures. Amorphous domains are seen: (a) to be characterized by a larger deviation of the rotational angle between thiophene rings out of coplanarity (±40° as opposed to ±20° for the crystal) and (b) to exhibit a significant percentage of anti-trans population

corresponds to a rotation of the ring by about 180° around the equilibrium crystal conformation, and was observed in all of our simulations with the amorphous P3HT systems irrespective of temperature. We also observed that the temperature dependence of the trans and anti-trans populations was the M


Macromolecules

Article

Although the conformation of the dihedrals is very similar in the pure crystalline and pure amorphous phases, deviations are recorded on the average size of the side hexyl branches. We see that at 300 K branches are more extended in the pure crystalline phase than in the semicrystalline one which is consistent with the larger gauche dihedral angle population in the latter case. More information about the microstructural changes taking place during the transition from the pure amorphous to the semicrystalline phase can be obtained by looking at the profiles of the radial distribution functions (RDFs) of different atom pairs in the simulated systems. Shown in Figure 20 are the partial RDFs for the intermolecular part of the ring center-ofmass (Figure 20a) and of the sulfur−sulfur pairs (Figure 20b) belonging to different chains, for the three characteristic temperatures: T = 300, 600, and 800 K. The two RDFs display the same behavior: with decreasing temperature there is a shift to smaller distances of the nearest neighbors and their intensity peak becomes more pronounced at 300 K (at around 3.97 Å). Moreover, as Figure 20b indicates, at 300 K new peaks start forming (the most significant at a short distance around 4.8 Å) reflecting a transition from a pure amorphous liquid-like behavior of the P3HT system at 800 and 600 K to a more structured glassy state where organization to longer characteristic distances is developing. Figure 20c compares the peaks at characteristic lengths of ordering for the RDF corresponding to the center-of-mass of the thiophene rings at 300 K for the glassy structure and its pure crystal analogue. The picture

same for the backbone thiophene ring dihedrals and the branch ones. A direct comparison of the branch size in the two phases at ambient temperature (T = 300 K) is displayed in Figure 19

Figure 19. Evolution of the instantaneous end-to-end distance of the hexyl branches for the pure crystalline and amorphous phases of the 20-3HT polymer at 300 K.

where the instantaneous square end-to-end branch distances for both crystalline and amorphous structures are plotted.

Figure 20. Radial distribution functions of (a) the center-of-mass of rings belonging to different chains (intermolecular part) and (b) between intermolecular sulfur pairs, for the amorphous 20-3HT system at T = 300, 600, and 800 K. (c) Comparison of the ring center-of-mass RDF between the pure crystalline and the pure amorphous phases of P3HT at 300 K. N


Macromolecules

Article

Table 3. MD Simulation Results for the Chain Packing Length, Square End-to-End Distance, and Persistence Lengtha lP (Å)

temperature (K)

a

± ± ± ±

crystal phase

370 335 300 225

11.3 11.0 10.8 10.3

0.4 0.4 0.5 0.4

amorphous phase

600 500 300

17.6 ± 0.8 16.8 ± 0.9 15.4 ± 0.8

⟨R2⟩ (Å2) 5574 5581 5586 5604

± ± ± ±

18 16 16 12

3204 ± 71 3187 ± 72 3171 ± 30

a (Å) 74.6 74.6 74.7 74.8

± ± ± ±

0.2 0.2 0.2 0.2

49.8 ± 2.4 49.6 ± 2.9 49.1 ± 2.5

The data are shown separately for chains in the crystalline and amorphous domains at several temperatures.

distribution of the CSCC dihedral (the torsion of the bond that connects two consecutive monomers along a P3HT chain) in crystalline and amorphous regions, respectively. The data for the simulation temperature of 300 K reveal that this bond holds significantly less all-trans population in the amorphous phase than in the crystalline one thus causing a more pronounced chain bending in the amorphous domains originating from fluctuations of this bond above and below the plane of the thiophene rings preceding and following it along the P3HT chain. This conclusion is further supported by inspecting the configuration of individual chains in the two phases. A typical example is shown in Figure 22 which reveals the greater flexibility of a chain in an amorphous zone at 300 K. Chain bending in the amorphous zones is very important in the light of a recent study by Crossland et al.55 on the efficiency of charge transport through intercrystalline amorphous domains, possibly through bridge chains. Our MD simulations show that chains in the amorphous phase are quite bent to significantly decrease the conjugation length. We have also analyzed the dynamics of the side hexyl segments for chains in the amorphous zone. The results are shown in Figure 23 and should be contrasted with the corresponding ones for hexyl side chain dynamics in the crystalline areas shown in Figure 10. Clearly, side chains undergo rather large-angle motions in the amorphous state. For chains in the corresponding crystalline domains, side chain motion was appreciable only for temperatures greater than 335 K. Again, by fitting the curves with KWW functions and computing the integral below the curves, we can get estimates of the corresponding correlation times; the results are reported in Table 4. By fitting the temperature dependence of these correlation times with a modified VFT function56 (Figure 23b)

provided by the RDF graphs is fully consistent with the visualization shown in Figure 16b for the structure of the amorphous phase of P3HT at 300 K. In an effort to get more insight into the structure and dynamics in the amorphous phase, especially in comparison with the corresponding information in the crystalline state, we have computed the packing length lp and the persistence length a for chains in the amorphous and crystalline domains and how their values change with temperature. The packing length can be computed from the MD predictions for the density of the system and the average square chain-end-to-end vector through: lp =

Mw ρNAv⟨R g 2⟩

(5)

where ρ denotes the density, Mw the molecular weight, NAv is the Avogadro number, and ⟨Rg2⟩ the average square radius of gyration of P3HT. The persistence length is computed by calculating the projection of the chain end-to-end vector R on the bond l1 connecting the S atoms in two consecutive odd (or even) P3HT monomers in the middle of the chain backbone, i.e., through: a=

⟨R·l1⟩ ⟨l1⟩

(6)

The simulation results are reported in Table 3 and show that chains are appreciably bent in the amorphous state, in complete contrast to chains in the crystalline domain which are mostly rod-like. To explain this, in Figure 21 we have plotted the

⎛ T ⎞ τc = τ0 exp⎜D ⎟ ⎝ T − T0 ⎠

(7)

(with D being a dimensionless parameter), we can compute the ideal glass transition or Vogel temperature for P3HT. This is the value of the parameter T0 in eq 7 and it came out to be between 245 and 285 K. The corresponding experimental estimate is 282.3 K.57

IV. CONCLUSIONS AND FUTURE PLANS In this work, we employed classical molecular dynamics simulations to examine the conformational behavior of RrP3HT molecules assembled into an ordered lamellar state of πorbital stacked conjugated backbones separated by regions of alkyl side chains. In particular, we focused separately on the pure crystalline and pure amorphous phases of Rr-P3HT as two different subsystem scenarios. Our purpose was to interpret the

Figure 21. Distribution of the CSCC dihedral for P3HT chains in the amorphous and crystalline domains at T = 300 K. O


Macromolecules

Article

Figure 22. Typical configurations of two randomly selected (a) amorphous and (b) crystalline P3HT chains from our MD simulations with the 203HT polymer at T = 300 K.

Figure 23. (a) Decay of the time autocorrelation function ⟨uBr(t)·uBr(0)⟩ for the end to end unit vector uBr of the 20-3HT hexyl branch in the amorphous phase, as a function of simulation temperature and (b) the temperature dependence of the correlation time τc for alkyl side group motion in the amorphous phase of 20-3HT. The dashed line denotes the best fit to the simulation data with a VFT function, eq 7 in the main text.

periodicity and therefore to higher charge carrier mobility. The π−π attractive interaction which is the main contribution for the P3HT self-assembly is not dramatically affected by temperature since the distance inside the crystal exhibits a minor fluctuation from 3.8 to 3.9 Å; therefore we expect temperature fluctuations not to radically affect the charge hoping in the π−π stacking direction. Furthermore, our calculations indicated that the planes containing the conjugated rings in P3HT are tilted by an angle about 20°, a result which is in accord with recent first principles density functional studies in combination with molecular dynamics simulations.29,30 As also stated in these publications, this rotation contributes substantially to the stability (the structure with tilt θ ≈ 20° corresponds to a local minimum of the potential energy) and affects substantially the electronic structure. As far as the dependence on the molecular weight is concerned we have come across a very significant result; we observed chain twisting for the longer crystals (30-3HT and 40-3HT) after a distance along the chain (around 8 nm) was exceeded, a result that is in perfect agreement with relative experimental findings.26 To the best of our knowledge no similar finding has been reported in simulation works thus far. In the case of the pure amorphous phase, well-relaxed configurations at the temperatures of interest were generated through a successful implementation of an extended cooling procedure cycle involving NVT-NPT MD relaxations. Our MD simulation results proved that the pure amorphous phase retains its melt configuration uniformity at high temperatures (∼600 K) but exhibits a clear transition to a glassy state at 300

Table 4. Values of the Characteristic Kohlrausch−Williams− Watts Relaxation Times, τKWW, Stretching Exponents, β, and Correlation Times, τc, of the KWW Functions Employed To Fit the Time Autocorrelation Function of the Hexyl Branch End to End Vector as a Function of Temperature

amorphous phase

temperature (K)

A

τKWW (ns)

800

1.0

∼0.41

0.22

∼26.3

600 500 300

1.0 1.0 1.0

∼4.40 ∼28.4 ∼1.4 × 105

0.23 0.24 0.19

∼165 ∼800 ∼29 × 106

β

τc (ns)

semicrystalline nature of Rr-P3HT by looking thoroughly into the structure of the two ideal subsystems composing it, leaving the study of the intermediate interfacial zone for a future investigation. In its pure crystalline form the system is characterized by a smectic-like structure with tilted side chains that do not interdigitate (no registry between layers) irrespective of temperature and molecular weight; this finding agrees with diffraction measurements in which the most stable polymorph, dictated as type I, was predicted to be noninterdigitated with tilted side chains.53 As the authors of ref 53 point out, the side chain attachment density of P3HT is too large to permit interdigitation compared with other new generation poly alkylthiophenes (like PBTTT and PQT) which display sufficient low attachment densities of their side chains allowing interdigitation thus leading to a three-dimensional crystalline P


Macromolecules

Article

points for subsequent higher-level simulations aiming to address the interplay between morphology (due to selfassembly) in the bulk of an organic monolayer or on electrodes and electron transport. By focusing first on small material samples, one can carry out quantum mechanics calculations60−67 to evaluate the key molecular parameters governing the rate of charge transfer between two adjacent molecules within a hopping regime, in amorphous and crystalline domains. Then, larger scale, kinetic Monte Carlo (kMC) algorithms can be utilized to describe charge transport through the entire material (containing amorphous and ordered domains) on the basis of the information for material morphology already obtained from the atomistic simulations and the precomputed charge hopping rates at the atomistic level for the system at hand.

K; at this temperature the system resembles more its semicrystalline nature rather than its pure amorphous phase. This picture was quantified by looking into the radial distribution function of the center-of-mass of rings belonging to different P3HT chains: at 300 K the structure of the RDF was found to deviate from that of a liquid-like melt system and started displaying peaks at higher characteristic distances. Currently, work is in progress to address significantly longer P3HT chain systems which will open the way to simulating semicrystalline P3HT samples at temperatures below the melting point. This is practically impossible to accomplish with MD methods which are known to be inefficient or even unfeasible because of the problem of long relaxation times, especially when highly ordered structures form at low enough temperatures. We mention for example that the MD simulation of a 48-chain 40-3HT system (P3HT chains made up of 40 thiophene rings) with the Dreiding force-field starting from the staggered configuration showed signs of chain twisting even at a low temperature (150 K) (see Figure 12) although the simulation was evolving too slow to let it go beyond 6 ns. In order to circumvent these MD drawbacks and be able to address longer P3HT chain systems, it seems imperative to resort to a nondynamic method such as Monte Carlo (MC). For systems like P3HT with large scale nanophase-separated structures at low temperatures, Metropolis Monte Carlo can tremendously help overcome all simulation difficulties associated with sluggish dynamics. We mention, in particular, MC algorithms based on a mixture of simple (such as reptation and generalized reptation, end-mer rotation, libration or flip, configurational bias, concerted rotation, and parallel rotation) and complex moves (such as end-bridging, double-bridging, double concerted rotation or d-CONROT, atom-identity exchange, and their variants to handle branch points)58 with which chain-like systems can be simulated in atomistic detail very efficiently. Thanks in particular to moves tailored for chain systems (such as the so-called family of chain-connectivity altering moves), equilibration in a MC process can be accelerated by a factor which can exceed the corresponding one obtained by MD by 4-to-5 orders of magnitude.59 We are currently redesigning all these moves for the P3HT chemistry, with the scope eventually to be able to simulate semicrystalline P3HT systems starting from randomly generated initial configurations. This will also open the way to extending our investigation to large, more realistic system sizes of P3HT similar to those employed in experimental studies or actually used in technological applications and to answer questions related with issues such as what is the dominant assembling mechanisms in P3HT, what is the resulting supramolecular arrangement of the formed two-dimensional foils or layers, what is the size of the interlamellar structure, and how regioregularity (stereoregular versus random-regular) affects the self-assembly and morphological properties of P3HT. On a longer time scale, such a powerful MC algorithm will allow us to study how differences in the chemical structure and/or architecture of the constituent chains can affect the morphological properties of a conjugated polymer. For example, we can carry out a comparative study of the morphological properties and local structure of poly-3hexylthiophene (P3HT) with poly[5,5′-bis(3-alkyl-2-thienyl)2,2′-bithiophene] (PQT) and with poly[2,5-bis(3-alkylthiophen-2-yl)thieno(3,2-b)thiophene] (PBTTT). The resulting MC configurations with the self-organized structures from the MC simulations will constitute ideal starting

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Telephone: +30-2610997398. Fax: +30-2610-965223. Author Contributions

The manuscript was written through contributions of both authors. Both authors have given approval to the final version of the manuscript. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS We acknowledge that the most part of the results in this paper has been achieved using the PRACE Research Infrastructure resource CURIE based in France at the Tre’s Grand Centre de Calcul (TGCC). We are also grateful to CINECA (Italy) and CSC-IT (Finland) for providing computational resources. This work has been funded by the MMM@HPC Project under the seventh Framework Programme of the European Commission within the Research Infrastructures (Grant Agreement Number RI-261594).

■

REFERENCES

(1) Salleo, A. Mater. Today 2007, 10, 38. (2) Heeger, A. J. J. Phys. Chem. B 2001, 105, 8475. (3) Virkar, A. A.; Mannsfeld, S.; Bao, Z. A.; Stingelin, N. Adv. Mater. 2010, 22, 3857. (4) Dimitrakopoulos, C. D.; Malenfant, P. R. L. Adv. Mater. 2002, 14, 99. (5) Faupel, F.; Dimitrakopoulos, C.; Kahn, A.; Woll, C. J. Mater. Res. 2004, 19, 1887. (6) Sariciftci, N. S.; Hoppe, H. J. Mater. Res. 2004, 19, 1924. (7) Kim, Y.; Cook, S.; Tuladhar, S. M.; Choulis, S. A.; Nelson, J.; Durrant, J. R.; Bradley, D. D. C.; Giles, M.; Mcculloch, I.; Ha, C. S.; Ree, M. Nat. Mater. 2006, 5, 197. (8) Yang, Y.; Li, G.; Shrotriya, V.; Huang, J. S.; Yao, Y.; Moriarty, T.; Emery, K. Nat. Mater. 2005, 4, 864. (9) De Boer, B.; Kroon, R.; Lenes, M.; Hummelen, J. C.; Blom, P. W. M. Polym. Rev. 2008, 48, 531. (10) Kanai, Y.; Wu, Z. G.; Grossman, J. C. J. Mater. Chem. 2010, 20, 1053. (11) von Wrochem, F.; Gao, D.; Scholz, F.; Nothofer, H. G.; Nelles, G.; Wessels, J. M. Nat. Nanotechnol. 2010, 5, 618. (12) Bao, Z.; Dodabalapur, A.; Lovinger, A. J. Appl. Phys. Lett. 1996, 69, 4108. (13) Raja, M.; Lloyd, G. C. R.; Sedghi, N.; Eccleston, W.; Di Lucrezia, R.; Higgins, S. J. J. Appl. Phys. 2002, 92, 1441. Q


Macromolecules

Article

(14) Mardalen, J.; Samuelsen, E. J.; Gautun, O. R.; Carlsen, P. H. Solid State Commun. 1991, 77, 337. (15) Prosa, T. J.; Winokur, M. J.; Moulton, J.; Smith, P.; Heeger, A. J. Macromolecules 1992, 25, 4364. (16) Ihn, K. J.; Moulton, J.; Smith, P. J. Polym. Sci. Polym. Phys. 1993, 31, 735. (17) Prosa, T. J.; Winokur, M. J.; McCullough, R. D. Macromolecules 1996, 29, 3654. (18) Hugger, S.; Thomann, R.; Heinzel, T.; Thurn-Albrecht, T. Colloid Polym. Sci. 2004, 282, 932. (19) Wu, Z. Y.; Petzold, A.; Henze, T.; Thurn-Albrecht, T.; Lohwasser, R. H.; Sommer, M.; Thelakkat, M. Macromolecules 2010, 43, 4646. (20) Malik, S.; Nandi, A. K. J. Polym. Sci., Polym. Phys. 2002, 40, 2073. (21) Joshi, S.; Pingel, P.; Grigorian, S.; Panzner, T.; Pietsch, U.; Neher, D.; Forster, M.; Scherf, U. Macromolecules 2009, 42, 4651. (22) Mccullough, R. D.; Tristramnagle, S.; Williams, S. P.; Lowe, R. D.; Jayaraman, M. J. Am. Chem. Soc. 1993, 115, 4910. (23) Dag, S.; Wang, L. W. J. Phys. Chem. B 2010, 114, 5997. (24) Pascui, O. F.; Lohwasser, R.; Sommer, M.; Thelakkat, M.; Thurn-Albrecht, T.; Saalwachter, K. Macromolecules 2010, 43, 9401. (25) Brinkmann, M.; Wittmann, J. C. Adv. Mater. 2006, 18, 860. (26) Brinkmann, M.; Rannou, P. Macromolecules 2009, 42, 1125. (27) Corish, J.; Feeley, D. E.; MortonBlake, D. A.; Beniere, F.; Marchetti, M. J. Phys. Chem. B 1997, 101, 10075. (28) Maillard, A.; Rochefort, A. Phys. Rev. B 2009, 79−115207. (29) Northrup, J. E. Phys. Rev. B 2007, 76−245202. (30) Lan, Y. K.; Huang, C. I. J. Phys. Chem. B 2008, 112, 14857. (31) Cheung, D. L.; McMahon, D. P.; Troisi, A. J. Phys. Chem. B 2009, 113, 9393. (32) Do, K.; Huang, D. M.; Faller, R.; Moule, A. J. Phys. Chem. Chem. Phys. 2010, 12, 14735. (33) Vukmirovic, N.; Wang, L. W. J. Phys. Chem. B 2011, 115, 1792. (34) Vukmirovic, N.; Wang, L. W. Nano Lett. 2009, 9, 3996. (35) Vukmirovic, N.; Wang, L. W. J. Phys. Chem. B 2009, 113, 409. (36) Ruhle, V.; Kirkpatrick, J.; Andrienko, D. J. Chem. Phys. 2010, 132−134103. (37) Darling, S. B. J. Phys. Chem. B 2008, 112, 8891. (38) Darling, S. B.; Sternberg, M. J. Phys. Chem. B 2009, 113, 6215. (39) Melis, C.; Colombo, L.; Mattoni, A. J. Phys. Chem. C 2011, 115, 576. (40) Daoulas, K. C.; Ruhle, V.; Kremer, K. J. Phys.-Condens. Matter 2012, 24−284121. (41) Mayo, S. L.; Olafson, B. D.; Goddard, W. A. J. Phys. Chem. 1990, 94, 8897. (42) MacKerell, A. D.; Bashford, D.; Bellott, M.; Dunbrack, R. L.; Evanseck, J. D.; Field, M. J.; Fischer, S.; Gao, J.; Guo, H.; Ha, S.; Joseph-McCarthy, D.; Kuchnir, L.; Kuczera, K.; Lau, F. T. K.; Mattos, C.; Michnick, S.; Ngo, T.; Nguyen, D. T.; Prodhom, B.; Reiher, W. E.; Roux, B.; Schlenkrich, M.; Smith, J. C.; Stote, R.; Straub, J.; Watanabe, M.; Wiorkiewicz-Kuczera, J.; Yin, D.; Karplus, M. J. Phys. Chem. B 1998, 102, 3586. (43) Gasteiger, J.; Marsili, M. Tetrahedron 1980, 36, 3219. (44) www.scienomics.com. (45) Marcon, V.; Raos, G. J. Phys. Chem. B 2004, 108, 18053. (46) Marcon, V.; Raos, G.; Allegra, G. Macromol. Theor. Simul. 2004, 13, 497. (47) http://lammps.sandia.gov/; LAMMPS Molecular Dynamics Simulator. (48) Roger W. Hockney, J. W. E. Computer simulation using particles; CRC Press: Boca Raton, FL, 1988; p 267. (49) Marchant, S.; Foot, P. J. S. Polymer 1997, 38, 1749. (50) Tashiro, K.; Kobayashi, M.; Kawai, T.; Yoshino, K. Polymer 1997, 38, 2867. (51) Yamamoto, T.; Komarudin, D.; Arai, M.; Lee, B. L.; Suganuma, H.; Asakawa, N.; Inoue, Y.; Kubota, K.; Sasaki, S.; Fukuda, T.; Matsuda, H. J. Am. Chem. Soc. 1998, 120, 2047. (52) Meille, S. V.; Romita, V.; Caronna, T.; Lovinger, A. J.; Catellani, M.; Belobrzeckaja, L. Macromolecules 1997, 30, 7898.

(53) Kline, R. J.; DeLongchamp, D. M.; Fischer, D. A.; Lin, E. K.; Richter, L. J.; Chabinyc, M. L.; Toney, M. F.; Heeney, M.; McCulloch, I. Macromolecules 2007, 40, 7960. (54) www.sigmaaldrich.com/chemistry.html; Sigma-Aldrich Catalog. (55) Crossland, E. J. W.; Tremel, K.; Fischer, F.; Rahimi, K.; Reiter, G.; Steiner, U.; Ludwigs, S. Adv. Mater. 2012, 24, 839. (56) Papadopoulos, P.; Peristeraki, D.; Floudas, G.; Koutalas, G.; Hadjichristidis, N. Macromolecules 2004, 37, 8116. (57) Ngo, T. T.; Nguyen, D. N.; Nguyen, V. T. Adv. Nat. Sci.: Nanosci. Nanotechnol. 2012, 3, 045001. (58) Alexiadis, O.; Daoulas, K. C.; Mavrantzas, V. G. J. Phys. Chem. B 2008, 112, 1198. (59) Mavrantzas, V. G.; Boone, T. D.; Zervopoulou, E.; Theodorou, D. N. Macromolecules 1999, 32, 5072. (60) Sirringhaus, H.; Brown, P. J.; Friend, R. H.; Nielsen, M. M.; Bechgaard, K.; Langeveld-Voss, B. M. W.; Spiering, A. J. H.; Janssen, R. A. J.; Meijer, E. W. Synth. Met. 2000, 111, 129. (61) Kline, R. J.; McGehee, M. D.; Kadnikova, E. N.; Liu, J. S.; Frechet, J. M. J. Adv. Mater. 2003, 15, 1519. (62) Vogl, P.; Campbell, D. K. Phys. Rev. Lett. 1989, 62, 2012. (63) Troisi, A.; Cheung, D. L.; Andrienko, D. Phys. Rev. Lett. 2009, 102−116602. (64) Troisi, A.; Orlandi, G. Phys. Rev. Lett. 2006, 96−086601. (65) Zade, S. S.; Bendikov, M. Chem.Eur. J. 2007, 13, 3688. (66) Grozema, F. C.; van Duijnen, P. T.; Berlin, Y. A.; Ratner, M. A.; Siebbeles, L. D. A. J Phys Chem B 2002, 106, 7791. (67) Baranovski, S. Charge transport in disordered solids with applications in electronics; John Wiley and Sons Ltd: Chichester, U.K., 2006.

R


All-Atom Molecular Dynamics Simulation of ... - ACS Publications

Recommend Documents