Multi-Scale Modeling of Bulk Heterojunctions for Organic Photovoltaic Applications

We discuss a variety of recent approaches to molecular modeling of bulk heterojunctions (BHJs) in organic photovoltaics (OPV). These include quantum chemical calculations of the electron donor and acceptor molecules (such as polythiophenes and fullerenes), molecular simulations of their interactions in atomistic detail, mapping between different levels of resolution, coarse–grained (CG) modeling of larger scale structure, as well as electrical modeling. These calculations give a holistic view of these systems from local interactions and structure up to morphology and phase behavior. We focus especially on the evolution of the mesoscale morphology that is characteristic of BHJs. The simulations are compared to experimental results, both structural and dynamic, wherever such results are available. It turns out that CG models, which are clearly necessary to reach the length scales for morphology development, can accurately capture the large scale structure of atomistic systems in wide areas of the mixture phase diagram. They can also be used to study the dynamic evolution of the microstructure of a BHJ in a system approaching the device scale. On the other hand, atomistic simulations can lead to a geometrically based understanding of local structure which is crucial for charge separation and transport. Donor/acceptor mixtures for use in OPV devices have been the subject of intense investigation for the past decade. Recent focus on polymer based solar cells has been accelerated by excitement about mixtures of polythiophenes, specifically poly(3–hexylthiophene) (P3HT) with fullerenes, specifically [6,6]–phenyl–C61–butyric acid methyl ester (PCBM) (Padinger et al., 2003; Sariciftci et al., 1992). The most interesting aspect of the P3HT/PCBM mixture is that the power conversion efficiency (PCE) of OPVs resulting from this mixture can vary by almost an order of magnitude (Ma et al., 2005b). This range of efficiency values comes from variability in the hole mobility in the polymer (Mihailetchi et al., 2006), crystallinity of the polymer (Erb et al., 2005), domain sizes (Hoppe & Sariciftci, 2006; Yang et al., 2005), and absorbance of the polymer in the visible range (Al-Ibrahim et al., 2005; Zen et al., 2004) as a Multi-Scale Modeling of Bulk Heterojunctions for Organic Photovoltaic Applications


Introduction
We discuss a variety of recent approaches to molecular modeling of bulk heterojunctions (BHJs) in organic photovoltaics (OPV).These include quantum chemical calculations of the electron donor and acceptor molecules (such as polythiophenes and fullerenes), molecular simulations of their interactions in atomistic detail, mapping between different levels of resolution, coarse-grained (CG) modeling of larger scale structure, as well as electrical modeling.These calculations give a holistic view of these systems from local interactions and structure up to morphology and phase behavior.We focus especially on the evolution of the mesoscale morphology that is characteristic of BHJs.The simulations are compared to experimental results, both structural and dynamic, wherever such results are available.It turns out that CG models, which are clearly necessary to reach the length scales for morphology development, can accurately capture the large scale structure of atomistic systems in wide areas of the mixture phase diagram.They can also be used to study the dynamic evolution of the microstructure of a BHJ in a system approaching the device scale.On the other hand, atomistic simulations can lead to a geometrically based understanding of local structure which is crucial for charge separation and transport.Donor/acceptor mixtures for use in OPV devices have been the subject of intense investigation for the past decade.Recent focus on polymer based solar cells has been accelerated by excitement about mixtures of polythiophenes, specifically poly(3-hexylthiophene) (P3HT) with fullerenes, specifically [6,6]-phenyl-C 61 -butyric acid methyl ester (PCBM) (Padinger et al., 2003;Sariciftci et al., 1992).The most interesting aspect of the P3HT/PCBM mixture is that the power conversion efficiency (PCE) of OPVs resulting from this mixture can vary by almost an order of magnitude (Ma et al., 2005b).This range of efficiency values comes from variability in the hole mobility in the polymer (Mihailetchi et al., 2006), crystallinity of the polymer (Erb et al., 2005), domain sizes (Hoppe & Sariciftci, 2006;Yang et al., 2005), and absorbance of the polymer in the visible range (Al-Ibrahim et al., 2005;Zen et al., 2004) as a result of changes in processing conditions.In essence, the efficiency of the device is controlled by the nanoscale morphology of the blend.Since these initial studies of P3HT/PCBM morphology, there have been numerous articles that show how morphology develops (Li et al., 2007) and how the morphology can be controlled (Moulé & Meerholz, 2009).A large variety of techniques for measuring aspects of morphology have been applied.X-ray crystallography techniques can observe crystal packing in semicrystalline samples, but amorphous samples are not measurable (Chabinyc, 2008).NMR can give local (nearest neighbor) information, but does not extend to large enough distances to resolve domain size (Nieuwendaal et al., 2010;Yang et al., 2006).Electron microscopy can image domain sizes, but suffers from poor signal-to-noise due to very little contrast between P3HT and PCBM (Andersson et al., 2009;Oosterhout et al., 2009;van Bavel et al., 2010;2008;Yang et al., 2005).Scanning probe techniques mostly give information about the top surface, though some charge transport information can be derived from specialized scanning probe techniques (Groves et al., 2010;Pingree et al., 2009;Reid et al., 2008).The vertical segregation of materials can be determined non-destructively using refraction techniques (including ellipsometry (Campoy-Quiles et al., 2008;Madsen et al., 2011), X-ray (Andersen, 2009;Germack et al., 2010), and neutron scattering (Lee et al., 2010;Mitchell et al., 2004)).Alternatively, the composition of exposed surfaces can be determined using X-ray photoelectron spectroscopy (Xu et al., 2009).The interior composition of a mixed donor/acceptor layer has also been determined by etching away material using an ion beam and then measuring dynamic secondary ion mass spectrometry (SIMS) (Björstrom et al., 2005;2007).The measurement tools that materials science has to offer have been applied in abundance to the determination of the device morphology.However, each measurement technique has sample preparation requirements that can make it difficult to answer desired scientific questions.For example, a transmission electron microscopy (TEM) sample must be removed from the substrate, which requires that the morphology does not change during that sample removal process.One problem that cannot be solved even with a combination of measurement techniques is the determination of nano-scale information about the molecular arrangement in amorphous domains.In fact, the only techniques that can provide complete molecular scale information about amorphous organic domains are computational.Computer simulations provide means to better understand the relationship between microscopic details and macroscopic properties in solar cells.Most importantly, simulations can be used to explain experimental observations and to elucidate mechanisms that cannot be captured experimentally.In this chapter, it is intended to provide an overview of relevant computational techniques along with advantages and limitations of these techniques with regard to OPV materials, and summarize the research that has been performed to provide insight into the molecular mechanisms in such materials.We also discuss the continuum drift-diffusion model and kinetic Monte-Carlo models used to describe the electric characteristics of OPV devices.

Modeling techniques
Computational modeling of soft matter in general and of OPV in particular is clearly not a "one size fits all" problem.Even though the computational expense of detailed modeling is continuously decreasing, it is still the limiting factor in addressing the size and time The local energy minimum of the structure is found by making small steps towards the local energy minimization.A global structural minimum can be obtained by varying the initial configuration and re-minimizing for each one until a statistically relevant number of starting points lead to a common deep structural minimum.The strength of MM calculations is that they can be used to obtain increased molecular detail for systems for which a partial crystal structure is available.For example, if the stacking distance between adjacent polymer chains is known, MM/EM calculations could provide details on the location of the side chains or the exact angle of backbones with respect to each other.EM calculations are meaningless without a starting configuration that is very close to the ideal because locally rather than globally minimized structures will be the result of the calculation.This means that non-crystaline structures cannot normally be determined using EM.Molecular dynamics (MD) calculations take kinetic energy (KE) into account and can provide structural insight that is not predetermined from the initial configuration.The movement of each particle is determined from the sum of the bonded and non-bonded forces from surrounding atoms.In many cases (especially on the coarse-grained level) also partially random forces are used to thermostat the system, leading to Brownian motion.Since KE is determined by setting a temperature, the structures formed are not the lowest PE, but an ensemble of all possible structures at a given temperature weighted to minimize free energy.MD are often performed from an ordered starting point that results from MM/EM calculations.MD generates (in the case of constant volume) a set of configurations representing a canonical ensemble in correct time order, i.e. dynamic information is available.The disadvantage of MD calculations is that equilibration (when the system only moves within the correct ensemble) may require a long time, so due to limited computation speed and resources, most systems must be kept small and one can only afford short equilibrations (≤ 100ns).MD is a step-wise process.At each time step, the forces on each atom are calculated.Then the Hamiltonian equations of motions (which in the typical case of Cartesian co-ordinates are mathematically the same as Newton's equations) are integrated over each time-step.In most simulations temperature and/or pressure is fixed.Thus, an algorithm called a thermostat and/or barostat is needed to correct temperature/pressure and then the step is repeated.This generates a time series of conformations representing the chosen statistical ensemble.Most molecular simulations are performed using atomistic or semi-atomistic models.The advantages are obvious: there is a unique and strong connection to the underlying chemistry and the models are in many cases easily available (can be downloaded from various websites).Atomistic models normally contain Lennard-Jones (LJ) interactions between non-bonded atoms in addition to a variety of bonded interactions.The most important bonded interaction is the dihedral as the dihedral energies are in the range of thermal energy whereas all other bonded interactions are much larger compared with kT.Quantum-chemistry is often used to obtain the dihedral interaction parameters.As proper dihedral interactions must be periodic with respect to 360 • , a Fourier series is typically used to describe the dihedral potential.Also, quantum-chemistry is often used to calculate partial charges, which are needed for electrostatic interactions.Most other interactions are fit to empirical data as no simple and unique theoretical technique is available.A large number of simulation packages (e.g.Gromacs, Amber, LAMMPS, DL_POLY, Tinker etc) have been developed for this kind of modeling and are either freely available or for only a nominal charge.For details the reader is referred to a number of books on simulations (Allen & Tildesley, 1987;Frenkel & Smit, 1996).

32
Third Generation Photovoltaics www.intechopen.com The macroscopic properties of semiconductors used in OPV active layers depend on a number of factors.It is important to relate the chemical structure of the components of an active layer to the overall behavior of the system.For OPV applications, it is of interest to study the structural development of active layer mixtures on time and length scales that are relevant to device studies in order to accurately aid in interpretation of experimental findings.However, the computational limitations of atomistic simulations only allow for systems of the size of a few nanometers to be studied.The length scale of interest is on the order of 5-10 nm for the purposes of studying charge transport, and on the order of 100 nm for analyzing the phase separation behavior of the active layer components.Coarse-graining is the process of systematically grouping atoms from the atomistic model into "super-atoms."This reduces the number of interactions and degrees of freedom in the system allowing for simulations to be analyzed over longer length and time scales.For the simulation of P3HT/PCBM, our group achieved a 100× increase in computational speed by using a coarse-grained (CG) model, which could be applied to photovoltaic device-scale systems (Huang et al., 2010).After this brief overview of the techniques we will now discuss a number of recent modeling approaches using atomistic, CG, and continuum models.

Atomistic scale modeling
Organic photovoltaics are characterized by charge carrier mobility, which in turn defines their device efficiency.In BHJ devices, the minimum requirement for high charge carrier mobility is the presence of percolating pathways which enable the holes and electrons to reach the electrodes producing a photocurrent.Thus, device performance is strongly dependent on system morphology, both at the molecular level and over larger length scales.There are a number of factors that can affect the morphology and their description hinges not only on atomistic, but also on coarse-grained, descriptions.Morphology is critically dependent on the underlying substrates, percolation, thickness of the thin film, and on the processing and annealing methods and conditions.The latter is a significant challenge to simulate, as basically in simulations we always try to approach equilibrium as close as possible, but often in OPVs, equilibrium is neither attainable nor actually desired.Atomistic level simulations play a critical role in elucidating detailed mechanisms at the nano-scale level.Such calculations eventually can lead to systematic improvement of performance of solar cells.Fullerene derivatives are important electron acceptors and have been studied with different models in varying degrees of detail (Andersson et al., 2008;Arif et al., 2007;Choudhury, 2006;Girifalco, 1992;Hagen et al., 1993;Kim & Tománek, 1994;Qiao et al., 2007;Wong-Ekkabut et al., 2008).Polythiophenes have only recently been studied by computer simulations (Akaike et al., 2010;Botiz & Darling, 2009;Cheung et al., 2009a;Curco & Aleman, 2007;Do et al., 2010;Gus'kova et al., 2009;Widge et al., 2007;2008).We will discuss below simulations of different organic semiconductor materials where we first discuss fullerenes and thiophenes and then discuss a few other commonly used materials.

Modeling of fullerenes as electron acceptors
Buckminister fullerene and its derivatives are the most widely established and effective electron acceptors in OPV.Their ability to accept electrons from commonly used donors, such as photoexcited conjugated polymers, on a picosecond time scale and their high charge mobility (Singh et al., 2007) make them appealing as electron acceptors.Since the first Fig. 1.Electron mobility for C 60 with different chain lengths at 300 K.The lines without the points superimposed had the site energy differences calculated from a first-order approximation to the site energy differences between the two hopping sites.The lines with points superimposed, are the calculated mobilities when the site energy differences have been forced to 0 meV (MacKenzie et al., 2010).Reprinted with permission from J. Chem.Phys.132, 064904, 2010.Copyright 2010 American Institute of Physics.demonstration of polymer/fullerene solar cells (Sariciftci et al., 1993), there has been a significant effort to improve their processability from solutions and to improve performance by optimizing their morphology and energy levels (Troshin et al., 2009).Molecular packing of the OPV active layer is greatly affected by the conformation of the side groups for both small molecules and polymers.MacKenzie et al. (MacKenzie et al., 2010) have studied the effect of charge mobility with different side chains on C 60 .MD has been performed to generate a realistic material morphology for a series of four C 60 derivatives.This series of four C 60 derivatives are formed by attaching a saturated hydrocarbon chain of lengths 0, 20, 40, and 80 carbons to the C 60 via a methano bridge.The addition of a side chain has been found to disrupt the optimal packing of C 60 .Furthermore, Figure 1 shows how increasing the hydrocarbon chain length reduces charge mobility.This may stem from the increased size of the functional group, which pushes the C 60 molecules away from each other and decreases the number of neighbors close enough to electronically interact.As a result, as the functional group size increases, the overall transfer rate decreases thereby reducing the charge mobility (MacKenzie et al., 2010).The effect of temperature on a thin film of C 60 has been studied by performing a Monte-Carlo calculation simulating the physical vapor deposition process (Kwiatkowski et al., 2009).The calculations performed at 298, 523, and 748 K reveal morphologies that seem quite different: crystalline regions are estimated to be around 2 nm, 4 nm, and 6 nm in length, respectively.This follows the same trend as the experimental measurements (Chen. et al., 2001;Cheng et al., 2003), even though the exact grain lengths are smaller than the measured values as it is computationally too expensive to achieve device scale equilibration.However, the calculated radial distribution function (RDF) for the three morphologies are quite similar to each other and to the single crystal RDF (Kwiatkowski et al., 2009).This can be attributed to the fact that C 60 is a rotationally symmetric molecule with isotropic intermolecular interactions (Girifalco, 1992).This indicates that these spheres can pack efficiently, that even in partially disordered films the molecules are well connected to each other (Kwiatkowski et al., 2009), giving rise to the observed high electron mobility of 8 × 10 −6 m 2 /(Vs) in thin films of evaporated C 60 (Haddon et al., 1995) .One of the main advantages of OPVs is the ease of fabrication, due to their solution processability.Recently, it has been shown that a highly soluble derivative of C 60 , phenyl-C 61 -butyric acid methyl ester (PCBM, Figure 2), performs better than C 60 in solution processed OPV devices (Yu et al., 1995).The addition of the functional side group renders the acceptor more soluble and increased the acceptor strength, thus leading to a higher open-circuit voltage (V OC ) in the resulting device (Hummelen et al., 1995).PCBM has a high electron mobility of 2 × 10 −3 cm 2 /(Vs) at room temperature (Mihailetchi et al., 2003) and since balanced hole and electron charge mobilities reduce space charge build-up and increase the filling factor (FF) of devices, donor polymers with higher hole mobilities are desired for use with PCBM.This functionalized fullerene is relatively inexpensive and it forms segregated phases with many common donors to form mixed layers with ideal morphology (i.e., a grain size separation on the order of the exciton diffusion length and a bicontinuous network formed from the two components) (Ballantyne et al., 2010;Ma et al., 2005b;Thompson & Frchet, 2008).In order to properly understand the elementary processes in photovoltaics, a correct description of the electron acceptor morphology is important.MD calculations (with the all-atom OPLS force field parameters) performed at 300 K by Cheung et al. (Cheung & Troisi, 2010) correctly predict the experimental low-temperature crystal structure (Rispens et al., 2003) of PCBM (Figure 3), in which the electron hopping is facilitated by the molecular arrangement.Figure 4 shows the RDFs between the fullerene centers of mass, the phenyl Fig. 3. PCBM system assembly at T = 300 K, showing only the fullerenes (top) and the side chains (bottom).Layering of PCBM molecules in a zigzag pattern is evident (Cheung & Troisi, 2010).Reprinted with permission from J. Phys.Chem. C 114, 20479-20488, 2010.Copyright 2010 American Chemical Society.ring centers of mass, and between the fullerenes and phenyl rings (Cheung & Troisi, 2010).The fullerene-fullerene and fullerene-phenyl RDFs at 300 K, 400 K and 500 K suggest at best a weak influence of temperature due to the rigid structure of the fullerene cage.For the fullerene-fullerene RDF (Figure 4(a)) the peak is around 10.2 Å, in agreement with experiments (Xu et al., 1993) and previous simulation work carried out by MacKenzie et al. (MacKenzie et al., 2010).The first shell in the phenyl-phenyl RDF is around 5 Å (Figure 4(b)).This is significantly larger than the π-stacking distance in aromatic organic crystals (r = 3.8 Å), which suggests that the aromatic moiety in the PCBM side chain does not play a major role in the charge transfer process (Cheung & Troisi, 2010).The tail of PCBM helps to further stabilize the crystal, by the formation of weak hydrogen bonds between aromatic rings and oxygen in addition to the van der Waals interactions (Ecija et al., 2007;Napoles-Duarte Fig. 4. RDFs for PCBM between, (a) Fullerene-fullerene (b) phenyl-phenyl (c) fullerene-phenyl, at T = 300 K (solid line, black), T = 400 K (dotted line, red), and T = 500 K (dashed line, green) (Cheung & Troisi, 2010).Reprinted with permission from J. Phys.Chem. C 114, 20479-20488, 2010. Copyright 2010American Chemical Society. et al., 2009).Since only the fullerene cage in PCBM is likely to be involved in charge transport, morphology related calculations may be used to describe the charge separation, transfer, and recombination processes in OPV BHJ.MD calculations have also been performed to examine the structural variation between PCBM molecules.Figure 5 shows the dihedral angle distribution for the phenyl group and ester group of the PCBM with respect to the fullerene.The probability distribution of φ 1 , the angle between the alkyl chain and phenyl ring, has two equally populated peaks at 90 • and 270 • (Figure 5(a)), indicating that the ring orients parallel to the fullerene surface to reduce steric hinderance.The side group alkyl chain predominately adopts a trans conformation (φ 2 = 180 • ).While at 300 K the population of the gauche conformations is around 15 %, this increases to around 20 % percent at 500 K (Cheung & Troisi, 2010).This is in agreement with crystallographic studies (Rispens et al., 2003).That is, this structure provides a good representation of the solid-state PCBM system.It is well-known that the PCE of a solar cell depends on V OC , which in turn is a function of the energy gap between the highest occupied orbital (HOMO) in the donor and the lowest unoccupied orbital (LUMO) in the acceptor.That is, increasing the LUMO in the acceptor improves the device PCE (Kooistra et al., 2007).An effective approach to increase the donor LUMO is to alter the fullerene cage by directly adding multiple substituents (Lenes et al., 2008).On the other hand, this increases the structural disorder of the fullerenes.A DFT study has been conducted (Frost et al., 2010) to study the balance between these two effects on device performance (Lenes et al., 2009).The authors of that work calculated the energy levels of all eight unique isomers of bis adduct of PCBM and ten representative isomers of tris-PCBM (Djojo et al., 1999).Indeed, a generally lower range of LUMO energies was calculated as a function of the isomers compared to mono-PCBM (Frost et al., 2010).Yet this large range of energetic disorder, in accordance with cyclic voltammetry measurements of the adducts and device current-voltage characteristics (Lenes et al., 2009), decreases the device FF imposing an adverse effect on performance (Frost et al., 2010).It is suggested that designing higher adduct fullerenes with shorter sidechains or hetero-adducts with distinct roles for each of the sidegroups or manipulating energy levels, will lead to fullerene derivatives with better performance as acceptors in solar cells without the need for expensive isolation of isomers (Frost et al., 2010).In that sense, atomistic level calculations can provide insight into the design and optimization of novel electron acceptors.

Atomistic structure and dynamics of polythiophene semiconductors
The electrooptical properties of semiconducting conjugated polymer films depend critically on molecular structure, doping, and morphology.Whereas the former two can be relatively well controlled during synthesis, the latter is far from well controlled or understood in experiments.Thus, accurate computer simulation models play an important role in elucidating the phase behavior of polymers.They can determine the most likely structures that a polymer will form and calculate the energetic price for assuming another structure.They aid the interpretation of experimental data, because particle positions can be tracked exactly during the course of a simulation.Poly(3-alkylthiophene)s [P3ATs] are an abundant class of polymer semiconductors.In particular, the hexyl derivative poly(3-hexylthiophene), P3HT (Figure 6), has become one of the most widely used electron donor materials in OPV.P3HT exhibits high charge mobilities (Sandberg et al., 2002).The charge mobility in polymer semiconductors is strongly dependent on the packing of polymer chains in P3HT films (Kline & McGehee, 2006;Li et al., 2005).Low molecular weight (LMW) (Meille et al., 1997) P3HT exhibits a preferred interdigitated structure, with the distance between the chains 3.8 Å perpendicular to the rings and 16 Å parallel to the rings (Kline et al., 2007;Prosa et al., 1996;Yamamoto et al., 1998).Radial distribution functions (RDFs), g(r), are a standard tool of quantitative structural characterization.RDFs describe how the atomic density varies as a function of distance from a reference atom.For a molecule containing several different atom types, partial RDFs can be defined between any two atom types.Cheung et al (Cheung et al., 2009b) have performed classical MD simulations at several different temperatures to study the microstructure of P3HT, using a force field previously employed in tetrathiophene calculations (Marcon & Raos, 2006;Marcon et al., 2006).This study on LMW P3HT (20 monomers per chain, molecular weight around 6700 Da) with an initially interdigitated configuration of three polymer layers depicts some important characteristics of P3HT.RDFs between several atoms in the same layer of polymer are illustrated in Figure 7.A common feature among all the RDFs is the peak at r = 4Å , which corresponds to the separation between monomers on the same chain.As temperature increases, it is expected that the interlayer lattice spacing will increase with the expansion of the alkyl chain.This is clearly seen in Figure 7 (c), where the position of the first peak grows with temperature.Such an increase in the interchain spacing with temperature is bound to decrease the P3HT charge mobility, because the probability for charge transfer between adjacent sites in the polymer matrix has an exponential dependence on distance.The sharp peaks at short distances in Figure 7  represents a less well-defined molecular structure (increased disorder) at high temperatures for LMW P3HT (Cheung et al., 2009b).Furthermore, at low temperatures and high temperatures, the system density shows a roughly linear decrease with temperature, with a distortion from linearity around 250 K between two linear regions.This is accompanied by a peak in the molar heat capacity, C p , and isothermal compressibility, k T , (Figure 8) around the same region.These features signify a conformational transition in the side chains, in accordance with X-ray measurements (Prosa et al., 1992): at low temperatures the ring-tail dihedral angle is around 100 • , but above Atomistic molecular dynamics simulations of P3HT and poly(2,5-bis(3tetradecyllthiophen-2-yl)thieno[3,2-b]thiophene) (PBTTT) at finite temperatures have been compared to investigate the nanoscale structural properties that lead to the higher measured hole mobility in PBTTT than in P3HT field-effect transistors (FET).Simulations of the two polymer melts show that the structural properties in PBTTT facilitate both intraand inter-chain charge transport compared with P3HT, due to greater degree of planarity, closer and more parallel stacking of the thiophene and thienothiophene rings, and possible interdigitation of the dodecyl side chains.X-ray diffraction studies have shown that PBTTT indeed forms interdigitated alkyl side chains (Brocorens et al., 2009;Kline et al., 2007).Thus, the crucial role played by the bulky dodecyl side chain and thienothiophene ring, respectively, in determining intra-chain and inter-chain structural order is clarified through these simulations (Do et al., 2010).In addition to determination of the interaction of a polymer with itself or with the fullerene, MD modeling can address anisotropic interactions such as the interaction of the polymer with a surface.For the most part, MD calculations of interactions between oligomers of polymers with a surface have been performed starting from a crystalline polymer structure that is then lowered onto the surface.MD simulations of thiophenes on surfaces have been recently reviewed (Gus'kova et al., 2009).Modelling of P3HT on a ZnO with the specific application of OPV devices was modelled by Saba et al. (2011) This group allowed a slab of P3HT 12-mers interact with a ZnO (10 10) surface at 2 K and determined that the structure formed by the P3HT is face on to the ZnO and also that the P3HT crystal spacing increases to lattice match to the ZnO.This article shows that surface properties can significantly influence the structure of polymers.The results, however, would be more believable if MD modeling had been performed at higher temperatures where the polymer had the opportunity to freely rearrange.It is clear that atomistic level calculations can accurately elucidate mechanistic details of the polymer matrix morphology.Due to the direct access to the microscopic scale that atomistic simulation affords, it can provide a vital link between molecular structure and charge mobility.Atomistic simulations run much more quickly for small molecules than for polymers due to the lower particle number and faster dynamics of the molecules near room temperature.There is a large body of literature involving atomistic modeling of other organic and organometalic compounds for evaporated photovoltaic devices.Discussion of this literature is beyond the scope of this book chapter.

Some other materials used in OPVs
Scanning force microscope (SFM) measurements on asymmetrically substituted poly[2-methoxy-5-(3',7'-dimethyloctyloxy)-1,4-phenylene vinylene] (OC 1 C 10 -PPV), also known as MDMO-PPV, shows spiraling chains and no aggregation.Both atomistic scale MD and Monte-Carlo calculations consistently portray this connected ring-forming structure formation of MDMO-PPV (Figure 9) (Kemerink et al., 2003).The bending force present in MDMO-PPV can be explained by an interaction, either attractive or repulsive, between the aliphatic side chains of successive monomer units.These seemingly connected rings could be due to the presence of a conformational or configurational defect.At such a defect, the orientation of the side chain is flipped from one side of the polymer chain to the other, and the bending direction is reversed (Kemerink et al., 2003).This molecular conformation and the weakly ordered stacking of the molecules, which result in a high disorder energy and poor π − π interaction, might be the origin of the poor transport properties in MDMO-PPV.The hole mobility of neat MDMO-PPV is around 5 × 10 −7 cm 2 /(Vs) (Blom et al., 1997;Martens et al., 1999;Vissenberg & Blom, 1999), whereas neat PCBM has an electron mobility that is 4000 times larger than that (Mihailetchi et al., 2003).Thus, charge transport in a BHJ OPV with an active layer of a MDMO-PPV and PCBM blend is expected to be strongly unbalanced.As the electrons do not neutralize the holes in the device, this results in an augmentation of space-charges and a limited photocurrent.However, hole mobility measurements in a blend of MDMO-PPV with an excess PCBM shows that hole mobility is increased.When mixed with PCBM, the ring formation in MDMO-PPV is hindered due to interactions between MDMO-PPV and PCBM (Figure 10).The change in morphology with an enhanced number of percolating pathways results in improved charge-transfer properties in the blend (Melzer et al., 2004).Hexabezocoronene (HBC) is a well-known discotic liquid crystal which self-organizes to form columns with strong π-orbital interactions between molecules within columns.Such materials, often referred to as an alternative to thin films of organic semiconductors, possess many characteristics that are important in solar cells (Funahashi, 2009): improved vertical mobility, self-annealing and self-assembling characteristics that give rise to larger domains thus reducing charge-trapping sites.Time-resolved microwave conductivity (TRMC) measurements have shown that hole mobility in HBC is well above 0.1 cm 2 V −1 s −1 and depends on side chains (van de Craats et al., 1999).Atomistic simulations (Andrienko et al., 2006;Kirkpatrick et al., 2008;Marcon et al., 2008) performed on several HBC derivatives with different types of side chains (Figure 11) illustrate the dependence of morphology on semiconductor side chains.Alkyl side chains are generally added to polymers to improve Fig. 9. Simulation of the surface morphology of MDMO-PPV using a Monte Carlo model.The inset illustrates the proposed side chain orientation.At a conformational defect (thick arrow), the sense of spiralling (thin arrows) reverses, leading to a morphology of connected rings (Kemerink et al., 2003).Reprinted with permission from Nano Lett. 3, 1191Lett. 3, -1196Lett. 3, , 2003. .Copyright 2003 American Chemical Society.(Melzer et al., 2004).Reprinted with permission from Adv. Func.Mater. 14, 865-870, 2004. Copyright 2004 Wiley.solubility and reduce the melting point.However, this increases the distance between separate polymer chains, resulting in lower charge mobility.As shown by both experiments and simulations, linear side chains show maximum packing order, while packing is significantly reduced for branched side chains.Consequently, long linear side chains show higher charge mobility than HBC with bulky side chains.Thus, whether it be a crystalline material or smectic liquid crystalline phase, the higher the degree of order, the higher the charge mobility of the material.

Coarse-grained modeling
A wide variety of different techniques for coarse-graining, especially for polymeric systems, have been developed (Baschnagel et al., 2000;Faller, 2004).The common main theme is to circumvent the problems associated with the small time and length scales in atomistic simulations by coarsening the model and decreasing the number of degrees of freedom.These mesoscale or CG models must conserve the chemical nature of the atomistic model.The interactions that describe the polymer properties can be determined for the coarse-grained system from a microscale reference simulation using, e.g., the Iterative Boltzmann Inversion (IBI) method (Reith et al., 2003).One groups atoms appropriately into "super-atoms" or CG beads using chemical intuition.The set of effective potentials to describe the structural distributions of the CG polymer chain can then be generated from a set of correlation functions obtained from a corresponding atomistic simulation.The IBI method is best described with an example involving a non-bonded potential, V 0 (r) (Faller, 2007;Reith et al., 2003).If we restrict ourselves to pair potentials non-bonded interactions can be fully described by a radial distribution function (RDF) which, describes how on average the particles in a system radially pack around each other.This radial packing illustrates the correlation between the packing of particles and the forces the particles exert on each other.Mathematically RDFs can be calculated by choosing a reference molecule in a system and a series of concentric spheres around it.RDFs are a measure of the number of sites in a sphere at distance r from the reference center divided by the number in an ideal gas at the same density: Here, g(r) is the RDF, n(r) is the average number of molecules in the shell which are counted based on the position of the center of mass of the molecule, d is the overall density, and 4r 2 ∆r is the volume of a spherical shell.The IBI aims to match the RDF from the coarse grained model onto the atomistic RDF by iteratively altering the CG potentials.In order to obtain an accurate reproduction of structural details in the CG model, it is important to determine the proper set of potentials.This equation defines the iterations, where, V i+1 (r) and V i (r) are the potential energies at the i + 1 th and i th iteration steps, respectively.The RDF of the coarse grained model at the i th step is described by g CG i (r) and is calibrated against the target atomistic RDF, g A (r).A reasonable initial guess for V 0 (r) is first required.The potential of mean force, F(r), based on the atomistic simulation will usually suffice: The iteration continues until the difference in the RDFs between the CG and the target is below a pre-defined tolerance.This method can be applied to any set of interactions by replacing the radial distribution function by the appropriate probability distribution and the potential by the correct interaction correlation.In dense systems, the interaction distributions are interdependent; thus, one cannot determine each potential separately.Instead, one normally performs the iteration on one of the potentials constant while keeping the rest constant.There must, in the end, be a readjustment of potentials after all are individually iterated.The speed at which this iterative process converges relies on the order in which one optimizes the set of potentials; it is better to start with potentials that are least affected by changes to the rest of the set (Huang et al., 2010;Sun & Faller, 2006;Sun et al., 2008).

Polymer/fullerene coarse-graining
Our group recently published the first systematic CG-MD model of OPV relevant compounds (Huang et al., 2010;2011), namely P3HT and the simplest fullerene, C 60 .C G simulations allow an increase in system size that could reasonably be simulated using a computer cluster to about 25 × 25 × 25 nm 3 (compared to about 6 × 6 × 6nm 3 atomistically).This larger volume allows the formation of morphological features approaching what is expected in a device.The larger volume also allows simulation of polymers with a molecular weight approaching typical experimental systems (9-18 kDa) (Brinkmann & Rannou, 2007;Schilinsky et al., 2005;Zen et al., 2004).In a follow up study, we simulated random mixtures of P3HT and C 60 that were equilibrated at high temperatures and cooled down to temperatures at which C 60 formed clusters.The coarse-grained model replicated many thermodynamic features that were physically expected.C 60 did not form large clusters in low MW P3HT but did so with higher MW P3HT.Also, C 60 did not form clusters until its concentration in P3HT reached a certain threshold (Fig. 12) (Huang et al., 2011).
For the first time a system approaching the size needed for domain formation was studied.The simulations clearly were able to demonstrate that fullerene forms disordered clusters as expected.But now the shape and size distribution of such clusters can be measured.It also turned out that the polymer conformations were much more heterogeneous than expected.In particular, the heterogeneity (e.g. in gyration radius and anisotropy) was found to vary with chain length.

Electrical modeling
Given a model for the BHJ morphology and its dependence on material properties and processing conditions, a theoretical understanding of the relationship between the morphology and charge-carrier and exciton transport is needed in order to predict the morphology dependence of device characteristics.A truly accurate theoretical description of charge-carrier and exciton transport must inherently be multi-scale in nature, because it must account for the anisotropic transport in conjugated polymers on the molecular scale (e.g. higher intra-versus inter-chain charge mobility), exciton dissociation and diffusion between electron donor and acceptor domains on the 10 nm scale, and charge-carrier transport within donor and acceptor domains on the 100 nm scale.

Continuum drift-diffusion models
The modeling of the electrical characteristics of OPVs has been heavily influenced by the more developed field of inorganic semiconductor physics.One common approach to device-scale electrical modeling of organic solar cells that has been translated almost directly from inorganic semiconductor field is the continuum drift-diffusion model (Sze & Ng, 2006).This approach involves solving the continuity equations for the electron and hole densities, n and p: assuming that the electron and hole current densities, J J J n and J J J p , consist of a drift term that is proportional to the gradient of the electrical potential ψ and a diffusion term proportional to the gradient of the carrier densities, Here e is the elementary charge, D is the charge-carrier generation rate (which is equal to the exciton dissociation rate for organic solar cells), R is the recombination rate, µ µ µ n and µ µ µ p 46 Third Generation Photovoltaics www.intechopen.comare the electron and hole mobilities respectively, and D n and D p are the electron and hole diffusion coefficients.Typically, the mobilities and diffusion coefficients are assumed to obey the Einstein relation, D n,p = k B Tµ µ µ n,p /e.In general, the mobilities and diffusion coefficients in anisotropic materials such as conjugated polymers are tensors, but they have been assumed to be scalars in most theoretical treatments.The drift-diffusion model is completed by the Poisson equation relating the electrostatic potential to the charge-carrier densities, where ǫ is the dielectric permittivity of the medium.
In contrast to inorganic semiconductor devices, in which light absorption produces free charge carriers directly, electrical modeling of organic photovoltaics must take into account the dynamics of (singlet) excitons and their finite probability of dissociating into free carriers.Some calculations have avoided explicit consideration of exciton dynamics by taking the charge-carrier generation rate to be a (constant) parameter that is chosen to fit device electrical characteristics (Maturová et al., 2009a;b).In the framework of the drift-diffusion model, the time evolution of the exciton density x is given by where µ µ µ x is the exciton mobility, G the exciton generation rate, and R d the exciton decay rate, which is usually assumed to have the form R d = x/τ x , where τ x is the exciton lifetime.The factor of 1/4 in front of the charge-carrier recombination rate in Equation 9accounts for the fact that only a fraction of charge-carrier pairs recombine to form singlet excitons (Buxton & Clarke, 2006;Shah & Ganesan, 2009).The recombination rate of free charge carriers is generally assumed to be of the bimolecular Langevin form, R = γnp, where γ = e µ p + µ n /ǫ, while the exciton dissociation rate D is often given by the Onsager theory of electrolytic dissociation (Onsager, 1934), generalized to account for electric field dependence by (Braun, 1984).More sophisticated approaches to the drift-diffusion model have also been employed, accounting for phenomena such as trapped charges (Hwang et al., 2009) and various forms of the density of states (MacKenzie et al., 2011).The exciton generation rate G is equal to the dissipation rate of optical power and is thus proportional to the absorbed light intensity.Some drift-diffusion calculations of BHJ solar cells have assumed an exponential attenuation of the light intensity in the active layer with attenuation rate proportional to the absorption coefficient (Buxton & Clarke, 2006;Shah & Ganesan, 2009), while others have used a more sophisticated optical transfer matrix approach (Moulé & Meerholz, 2007;Pettersson et al., 1999) to account for interference effects between ingoing and outcoming light waves (Kotlarski et al., 2008;Nam et al., 2010).
Interference results in significant variations in the light intensity in thin-film devices for layer thicknesses on the order of the wavelength of light.Some calculations have simply assumed a constant generation rate within the active layer (Koster et al., 2005;Maturová et al., 2009a;b).While this may seem a drastic approximation, assuming a constant generation rate or using the optical transfer matrix method with the same total generation rate has been shown to make little difference to the the electrical characteristics of MDMO-PPV:PCBM BHJ solar cells (Kotlarski et al., 2008).
To solve the drift-diffusion equations, boundary conditions at the electrode interfaces for the electrical potential and charge-carrier densities (or current densities) must be specified, which in general depend on the applied voltage, electrode work functions, and electron donor and acceptor HOMO and LUMO levels (Koster et al., 2005;Shah & Ganesan, 2009).For charge transfer at internal interfaces between donor and acceptor domains, a transfer rate that depends exponentially on the energy difference between the donor and acceptor HOMO levels (for hole transport) or LUMO levels (for electron transport) is often assumed (Ruhstaller et al., 2001).Device characteristics are generally calculated under steady-state conditions, in which the time derivatives on the left-hand sides of Equations 4, 5, and 9 are set to zero.The parameters used in the drift-diffusion model, such as the charge-carrier mobilities, exciton lifetime, dielectric permittivity, electrode work functions, and donor and acceptor HOMO and LUMO levels, are generally taken from experimental measurements of the pure donor, acceptor, and electrode materials.These properties are therefore assumed to be isotropic and spatially invariant, except that differences between the values in donor and acceptor domains are sometimes taken into account.In some cases spatial variation of the carrier mobility is introduced in the form of a Poole-Frenkel electric field dependence, µ n,p = µ 0n,p exp E/E 0n,p , where E = |∇ψ| is the magnitude of the electric field and µ 0n,p and E 0n,p are constants (Buxton & Clarke, 2006).In the simplest examples of drift-diffusion modeling of BHJ solar cells, the structure of the BHJ has been ignored entirely and the active layer has been assumed to be completely uniform (Hwang et al., 2009;Koster et al., 2005;Sievers et al., 2006).In this case, the model simplifies considerably, with the electric field and charge carrier densities depending only on one spatial dimension.Such models can be useful for understanding the dependence of device characteristics on parameters such as the active layer thickness, injection barriers, and carrier mobilities.However, they completely ignore the sensitivity of device characteristics to the BHJ morphology, which is demonstrated by the strong dependence of solar cell efficiency on processing conditions for the same donor/acceptor blend ratios (Saunders & Turner, 2008).The BHJ morphology has been accounted for in some drift-diffusion calculations (Figure 13) by assuming a regular two-dimensional array of interdigitated donor and acceptor domains characterized by a domain width or widths (Maturová et al., 2009a;b;Nam et al., 2010).The material properties within the domains in these calculations have been taken to be uniform.These models allow the dependence of the device characteristics on the level of phase separation in the BHJ to be determined but they assume a BHJ morphology rather than determining it from some structural model.Within the framework of drift-diffusion models, some calculations have accounted for the effect of material properties and thermodynamic conditions on the heterojunction morphology formed and the consequences for electrical characteristics.Buxton et.al. (Buxton & Clarke, 2006) used a Flory-Huggins Cahn-Hilliard model to determine the morphology of an active layer consisting of a donor-acceptor block copolymer under various conditions and a two-dimensional drift-diffusion model as described above to calculate device properties.Adding another level of complexity (Shah & Ganesan, 2009), used self-consistent field theory (SCFT) to calculate the morphology of a donor-acceptor rod-coil block copolymer, accounting for the orientational anisotropy of the rods and the resulting charge mobility anisotropy.With this model, they found a significant interplay of domain size and anisotropy in optimizing device characteristics.

48
Third Generation Photovoltaics www.intechopen.comFig. 13.Low field current-voltage characteristics of BHJ solar cells with different length scales of phase separation.Letters and numbers in the legend refer to solvent [toluene (t) or chlorobenzene (cb)] and spin-casting speed (in rpm), respectively.The symbols are the experimental data, and the lines were calculated using the parameters employed in the numerical modeling.The inset shows the sample layout for the numerical model.Each slab has a thickness d and repeats with a dimension L where L d and L a represent the width of the donor and acceptor slabs, respectively (Maturová et al., 2009b).Reprinted with permission from Nano Lett.9, 3032-3037, 2009.Copyright 2009 American Chemical Society.
The structural models of Buxton & Clarke (2006) and Shah & Ganesan (2009) assumed generic material properties rather than considering specific realistic systems.They were also restricted to BHJs consisting of block copolymers, for which the equilibrium morphologies, which are readily obtained by methods such as SCFT, consist of phase-separated donor and acceptor domains.For typical BHJs consisting of donor and acceptor materials that are not chemically bonded and for which the equilibrium morphology is two completely phase-separated donor and acceptor domains (i.e. a bilayer), such approaches are less suitable.In these cases, an accurate means for calculating non-equilibrium morphologies is required.Non-equilibrium continuum field theories of polymer mixtures or polymer/fullerene mixtures are not well developed, although progress has been made on this front (Ceniceros et al., 2009).Molecular dynamics simulations, such as those described above for atomistic and coarse-grained models, provide a well-established means for obtaining non-equilibrium morphologies.In particular, coarse-grained simulations like those described in the previous sections provide a means of obtaining non-equilibrium BHJ morphologies on length scales close to the device scale.Continuum drift-diffusion models of charge-carrier and exciton cannot, however, be directly combined with BHJ morphologies from molecular simulations to study the morphology dependence of device characteristics.The next section describes how molecular dynamics models can be used to describe charge-carrier and exciton transport.

Kinetic Monte Carlo models
At the other extreme to continuum drift-diffusion modeling, charge-carrier transport has been studied on the atomic scale using quantum chemical calculations and time-dependent perturbation theory (Brédas et al., 2004;Cheung & Troisi, 2008;Coropceanu et al., 2007).Such calculations are computationally intensive and are generally restricted to systems on the nanometer scale and therefore cannot be used to study transport between donor-acceptor domains on the 10 nm scale, let alone device-scale transport.Vukmirovic et.al. (Vukmirović& Wang, 2008;2009) have developed a multi-scale ab initio method of simulating charge-carrier transport that provides a systematic means of combining density functional calculations of small-scale local structural motifs to obtain the electronic structure and charge-transport properties of organic semiconductors up to the 100 nm length scale.But this method has only been applied to studying pure materials so far.The equations of time-dependent perturbation theory can be simplified in the case of weak electronic coupling between donor and acceptor sites.
In the limit of weak electron-phonon coupling and low temperature, the charge-transfer rate can be described by the Miller-Abrahams formalism (Coropceanu et al., 2007), in which the hopping rate from site i to j is where ν is the attempt-to-escape frequency, a 0 the localization radius, r ij is the hopping distance, and E i and E j are the energies of sites i and j respectively.In the limit of strong electron-phonon coupling and high temperature, the charge-transfer rate can be described by Marcus theory (Coropceanu et al., 2007;Nelson et al., 2009), where J is the electronic coupling between initial and final states (which depends, among other things, on the hopping distance r ij ), λ is the reorganisation energy, and ∆G is the difference in free energy between the initial and final states, which is often assumed to be equal to the difference in energy (E j − E i ).In the weak-coupling limit, exciton transport can be described by Förster resonant energy transfer (FRET) between sites, which, in its simplest form, gives an exciton transfer rate that is inversely proportional to the sixth power of the hopping distance r ij and decays exponentially with the energy difference (E j − E i ) between the sites (Watkins et al., 2005).The parameters in the transport equations can be estimated from quantum chemical calculations or experimental measurements (Coropceanu et al., 2007;Nelson et al., 2009;Westenhoff et al., 2006;2005).Such models of exciton and charge-carrier hopping between discrete sites have been used to study exciton and charge-carrier dynamics in BHJs, which can be simulated with a kinetic Monte Carlo algorithm (Frost et al., 2006;Groves et al., 2009;Meng et al., 2010;2011;Watkins et al., 2005).These models can account for the molecular nature of the active layer, the morphology of the BHJ, and the anisotropy of exciton and charge transport in a straightforward fashion.However, all implementations of such models so far (Frost et al., 2006;Groves et al., 2009;Meng et al., 2010;2011;Watkins et al., 2005) have treated the 50 Third Generation Photovoltaics www.intechopen.comdependence of the BHJ morphology on material properties and thermodynamic conditions in an approximate, albeit physically motivated, fashion: electron donor and acceptor sites were assumed to occupy sites on a cubic lattice and their interactions tuned to produce varying levels of phase separation after donor and acceptor sites were moved according to a Monte Carlo algorithm.Groves et.al used a kinetic Monte Carlo model in a cubic lattice to coarsen a random mixture of donor and acceptor sites into extended domains by allowing high energy interfacial sites to "exchange" position with lower energy sites.Once the coarsened morphology was generated, a drift-diffusion model was used to electrically model the device.When energetic disorder of isolated or smaller domains/sites was used to reduce charge mobility locally, the effects of local morphology on electrical properties could be modelled (Groves et al., 2009).Charge mobility was shown to be the most important determinant of efficiency using this model.
Recently, the Marcus charge-hopping model has been coupled with morphologies obtained from molecular dynamics or Monte Carlo simulations of more realistic molecular models to study charge transport in various physical pure materials, including hexabezocoronene (Kirkpatrick et al., 2007) and fullerenes (Kwiatkowski et al., 2009;MacKenzie et al., 2010;Nelson et al., 2009).But so far, such hopping models have not been combined with realistic structural models of electron donor/acceptor mixtures to study exciton and charge-carrier transport in BHJs.The atomistic and coarse-grained structural models described in the previous sections are amenable to such a treatment.By back-mapping the coarse-grained simulation configurations on to an atomistic description (Baschnagel et al., 2000) and estimating exciton and charge transport parameters from quantum chemical calculations of small atomistic systems as functions of relevant structural parameters, molecular-scale exciton and charge-carrier transport over an entire coarse-grained system can be studied.Such a strategy could be used to develop a multi-scale model of exciton and charge-carrier transport that accurately accounts for the BHJ morphology and its relationship to electrical characteristics from the molecular level up to the device scale.

Conclusions
It is obvious that the realistic description of OPVs is a true multi-scale problem.Therefore, the question of which technique is the best to describe the system is meaningless.We need all of them and much work needs to be done to combine them in a meaningful manner.We will always need to identify the right time and length scale for any given problem and then select the model to describe it, never the other way round.It may be tempting to take a model that has been validated in other contexts and just apply it.But we always have to make sure that it contains the right physics and chemistry and that it can realistically describe the relevant aspects of the system.This chapter has clearly been constrained to a small number of problems.But it should have given the reader a flavor of what is possible and how to approach modeling organic photovoltaics in a meaningful manner.

Fig. 12 .
Fig. 12.(a) Snapshots of configurations from atomistic and coarse-grained simulations of small systems with P3HT:C 60 = 1.85:1 w/w with N mono = 12 at 550 K and 1 atm(Huang et al., 2010).A single molecule of each type is highlighted.(b) C 60 aggregation in a simulation containing P3HT:C 60 = 1.27:1 w/w with N mono = 48.The largest C 60 cluster is highlighted in blue and all other particles in the system are shown as dots(Huang et al., 2011).Reprinted with permission from J. Chem.Theor.Comp.6, 527-537, 2010 and Fluid Phase Equilib.302,  21-25, 2011.Copyright 2010 American Chemical Society and 2011 Elsevier, respectively.