in crystalline PCBM-like fullerene derivatives : a comparative computational study †

Dipartimento di Chimica, Materiali e Ingeg Milano, via Mancinelli 7, IT-20131 Milano, +39-02-23993180; Tel: +39-02-23993051 Max-Planck-Institut für Kohlenforschung, K an der Ruhr, Germany † Electronic supplementary information mathematical derivation of the free e numerical data on the Huang–Rhys f Marcus–Levich–Jortner rate expressions, dependence of the mobility. An archive each crystal structure allowing 3D-visuali of the hopping probabilities derived (F3-network). See DOI: 10.1039/c4tc00502 ‡ Present address: Laboratory for Chem Mons, Belgium. Cite this: J. Mater. Chem. C, 2014, 2, 7313


Introduction
The synthesis of [6,6]phenyl-C 61 -butyric acid methyl ester (PCBM) was rst reported by Wudl and coworkers in 1995. 1 Today, aer almost two decades, PCBM is still the most popular n-type (i.e.3][4][5] It is likely that this enduring success would have surprised even PCBM's creators since, according to their original statements, 1 at that time they were simply aiming at a soluble fullerene derivative which could be produced easily and cheaply, and be readily adapted for a variety of purposes through further functionalization reactions.Several attempts have been made to nd good alternatives to PCBM, [6][7][8] sometimes using relatively sophisticated "molecular design" approaches (including ab initio calculations, for example, to select candidates with suitable electron affinities), but none of them has really taken hold.As a result, it is still very pertinent to ask: what makes PCBM special?Liu and Troisi have recently come up with one possible answer to this question: 9 unlike most other n-type materials, fullerene derivatives such as PCBM have anions with low lying excited states which may take part in charge separation at the interface with the donor polymer (exciton dissociation).This is interesting and useful, as it readily suggests a wide array of potential replacements for this material.On the other hand, this or any other design rule based on single-molecule properties cannot account for the full story, as they neglect all aspects related to molecular interactions and solid state organization, which are obviously crucial for properties such as processing and long-term stability, as well as charge and energy transport. 10everal experimental studies have evidenced that the structural organization of PCBM in photoactive blends ranges from largely amorphous, to nanocrystalline, to highly crystalline structures, [11][12][13][14] in which the solvent is thought to be absent. 15he relationship between the structural organization of PCBM molecules and their charge transport properties (electron mobilities, in this case) clearly represents a key issue.7][18] Indeed, in a previous study, 19 we addressed the rst part of this program by MD simulations of the structural organization of PCBM.Starting from the only two experimental crystal structures which were available at the time, 20 respectively for a 1 : 1 co-crystal with ortho-dichlorobenzene (DCB) and a 2 : 1 co-crystal with monochlorobenzene (MCB), we developed several models for amorphous, solventfree PCBM by means of different combinations of solvent abstraction, heating and cooling.On the basis of these simulations, we concluded that a key property of PCBM is that, even in the amorphous state, it can form well-developed threedimensional networks of electronically coupled fullerene molecules.Later on, Tummala et al. 21extended this study by carrying out further MD simulations and computing electron mobilities by kinetic Monte Carlo (KMC) simulations.More recently, other crystal structures have been added to the original pool, including the rst solvent-free PCBM structure reported so far 22,23 (these two papers discuss the same crystal structure, obtained respectively from powder and single-crystal X-ray diffraction) and that of the closely related ThCBM (1-thienyl-C 61 -butyric acid methyl ester), 24 which is obtained from PCBM by replacing the phenyl with a thienyl group (see Fig. 1).Our group has identied also other PCBM polymorphs, including one containing dichloromethane (DCM). 25In parallel with these structural investigations, other authors have performed computational studies of the transport properties of fullerenebased materials.MacKenzie et al. 26 combined MD and KMC simulations, demonstrating that longer side chains tend to interfere with electron transport.More recently, some authors have cast doubt on the applicability of the hopping models, arguing that the charge-localized picture which underlies this model is in fact inadequate for PCBM.Cheung and Troisi 27 highlighted an unusual distribution of thermally accessible localized and delocalized states, which cannot be mapped onto standard models of transport in disordered media.Blumberger and coworkers 28 computed electron transfer (ET) rates using a semi-classical rate expression including an adiabatic correction, and pointed out that the activation free energy for hopping can vanish for some strongly coupled fullerenes.
The general aim of the present work is to contribute to the understanding of the relationship between the solid-state organization and the charge transport properties of fullerene derivatives.Are electron mobilities within these fullerene-based materials highly sensitive to the formation of a specic polymorph?Is charge transport more anisotropic in some forms than in others?Could this be hampered by the possible inclusion of solvent molecules?Do the low-lying virtual orbitals characterizing these materials play a role in charge transport?As our starting point, we have taken all the presently available experimental crystal structures of PCBM and ThCBM, which may or may not include solvent molecules.The latter is not as widely used as PCBM, but it is nonetheless very important from our point of view because of its high electron mobility (2 cm 2 V À1 s À1 along one direction), as measured by time-resolved microwave conductivity on the same single crystals used for its structural determination. 24Electron mobilities have been evaluated using different versions of the hopping model for charge transport, despite its known limitations for fullerene-based materials. 27,28Compatible with these general aims and approach, we have tried to be as thorough as possible by checking the effect of different variables and computational assumptions.In particular, we have investigated: (1) the level of theory in the calculation of the electronic couplings (density functionals with different basis sets and semiempirical methods), (2) the intramolecular reorganization energies of different derivatives, (3) the effect of surrounding molecules on the site energies, (4) the adiabaticity of the electron hopping events, and (5) the complexity of the percolation network entering the KMC simulations, including the contribution of the fullerenes' low-lying virtual orbitals to ET.

Computational approach
We have adopted the approach reported in ref. 28 for the evaluation of the ET rates.On the basis of localized electronic states, we computed ET rates within different crystal structures using the same semi-classical rate expression and veried the non-adiabatic character of the different ET directions by calculating activation free energies.Note that the crystal structures were "frozen", with all atoms at their experimental (average) positions.The inclusion of thermal uctuation effects would have required extensive MD simulations of each crystal polymorph.The so-calculated ET rates were then used to perform KMC simulations to derive electron mobilities.For these simulations, percolation networks of various complexities were considered, including higher energy unoccupied orbitals on the fullerene derivatives.

Electron transfer rates
The ET rates were calculated 29 using the following semi-classical rate expression based on the Landau-Zener (LZ) treatment of non-adiabatic transitions 28,30 (see the ESI † for a mathematical derivation): with and In this expression, k el is the electronic transmission coefficient, n eff is the effective nuclear frequency, DE ‡ na is the nonadiabatic energy barrier and D is a correction factor relating DE ‡ na to the adiabatic energy barrier are respectively the electronic coupling and site energy difference between the initial and nal electronic states involved in the ET reaction, and l is the total reorganization energy.b is 1/k B T, where k B is the Boltzmann constant and T is the absolute temperature.
We shall see that under specic conditions these equations predict negative adiabatic "barriers".In such cases, we have still applied eqn (1) even though it can lead to ET which is faster than the vibrational time scale.This is an unphysical outcome of the localized hopping description, which however cannot be easily avoided at present.
The LZ approach was chosen to check the adiabaticity of ET directions, via the calculation of adiabatic barriers, consistently with other theoretical studies on PCBM. 28However, as many other workers in the eld adopt the Marcus (non-adiabatic limit of the LZ expression) or the Marcus-Levich-Jortner (MLJ, quantum-mechanical treatment of non-adiabatic transfer reactions including tunneling effects) expressions, we also computed ET rates using them.In order to keep our discussion mainly focused on the comparison of the different crystal structures, the explicit expressions and the numerical results obtained with these rates are given in the ESI.† Overall, we can say that the electron mobilities deriving from the different expressions are comparable in magnitude and they tend to be in the order m Marcus < m LZ < m MLJ .The reasons for this are also discussed in the ESI.† Thus, even though it would be hard to say whether the LZ approach is really superior to the other two, it has at least the additional advantage of interpolating between them.

Electron transport simulations
KMC simulations were performed adopting the Bortz-Kalos-Lebowitz (BKL) 31 algorithm to propagate a single charge carrier within the different crystal structures.For this purpose, we have developed an R-package 32 for KMC simulations of charge transport, which has been released on "Comprehensive R Archive Network" (CRAN). 29The package also computes the necessary rate constant (LZ, Marcus or MLJ) from their basic quantum chemical ingredients.The mobility tensor m for each structure was calculated from the expression ṽ ¼ À m F, in which F and ṽ are the applied electric eld (| F| ¼ 10 6 V cm À1 ) and dri velocity vectors.The minus sign accounts for the fact that the charge carriers are negatively charged.Our choice for the electric eld strength is discussed in the ESI, † which contains some numerical results for the eld-dependence of the mobility in one of the crystal structures.Independent simulations with the electric eld along the positive and negative x-, y-and z-directions were performed to build the mobility tensors.Each simulation lasted 10 7 hopping events and was repeated sixteen times with different random numbers (implying different initial positions of the charge).All our simulations were conducted by assuming that hopping events are restricted to neighbouring fullerene molecules, neglecting hopping to/from solvent molecules since these were found to display very high reorganization and site energies.Average hopping probabilities between each pair of molecules were extracted from these simulations and visualized.
The LUMO (L0), LUMO + 1 (L1) and LUMO + 2 (L2) orbitals are degenerate in C 60 and they are relatively close in energy also in its derivatives PCBM and ThCBM.For both of them, DE L0L1 z 0.06 eV and DE L1L2 z 0.22 eV, while DE L2L3 z 0.9 eV (B3LYP/6-31G* calculations).Since the thermal energy is k B T ¼ 0.026 eV at room temperature, in principle they could also participate in the electron transport.Somewhat loosely, we may speak of "excited state conduction", even though L1 and L2 are just oneelectron functions which do not correspond to any true excitation of the molecules or their anions.In our view, their true role is to provide some additional exibility in the description of the electrons' motion.We shall see that this partly compensates for the lack of thermal uctuations in our structural models, by allowing some hopping events also among molecular pairs which would be forbidden on the basis of the L0 orbitals.
In our KMC simulations, for each neighbouring pair of molecules (dened on the basis of the coordination numbers reported in Table 1), an electron can be transferred along nine different pathways, related to the interaction between different orbitals centered on each fullerene (L0L0, L0L1, L0L2, L1L0, L1L1, L1L2, L2L0, L2L1 and L2L2).The results from this "F3network" (where F indicates the full set of fullerenes, and 3 the number of orbitals per molecule) represent our standard for comparison with other, more approximate mobility simulations.First and foremost, they can be compared with more conventional KMC simulations based on a smaller network containing only the L0 orbitals of all molecules ("F1-network").Within the lower-symmetry triclinic structures (see below), nonequivalent molecular sites can display signicantly different L0, L1 and L2 energies.To determine to what extent the higher energy molecules contribute to ET, two further reduced percolation networks have been considered in the KMC simulations.The R3-network contains only the lower energy molecules, keeping all their three orbitals.The R1-network contains the same reduced set of molecules, keeping only their L0 orbital.

Electronic couplings and site energy differences
The electronic couplings and site energies were calculated by the fragment orbital approach described in ref. 33 and 34, using the ORCA 35 and cclib 36 packages.In this approach, for a single calculation the entire system is reduced to a dimer and the charge localization is mimicked by using the frontier orbitals of two isolated molecular fragments.In the case of an electron (hole) transfer, the Hamiltonian of the system is built in the basis of the LUMOs (HOMOs) of the monomers constituting the dimer.Aer re-orthogonalization, 37 the site energies and transfer integral can be extracted respectively from the diagonal and off-diagonal terms of the transformed Hamiltonian.Thus, we extracted dimers for each ET direction within the different crystal structures and, using neutral state fragment-orbitals, electronic couplings and site energies were calculated at the DFT or semiempirical levels.We point out that, according to a recent paper with high-level hole transfer calculations on a series of small dimers, 38 such fragment-based approaches can be expected to underestimate the electronic couplings by 20-30%.
In the approach described above, site energies are calculated on molecular pairs and, since these depend on the dimer used for the calculation, they cannot be considered as true site properties.Furthermore, a calculation based only on dimers will tend to underestimate polarization effects. 39Another possible route may consist of using micro-electrostatic (ME) methods, as recently done to investigate energetics at the P3HT/ PCBM interfaces 40 and at the interface between donor-acceptor discotic liquid crystals. 41However, for the present work, ME methods cannot be used as site energies have to be computed for different virtual orbitals.Thus, to obtain better dened and more realistic values, we considered larger clusters centered on each site.To calculate the L0, L1 and L2 energies on a given molecule, we generalized the fragment orbital procedure by building a Hamiltonian for a cluster made up of that molecule with its L0, L1 and L2 orbitals, and its rst coordination sphere with its L0 orbital.These site energies will be compared to those from pair calculations in the Results and discussion section.

Reorganization energies
The total reorganization energy involved in the electron transfer event investigated here is composed of an intra-molecular contribution l i and an outer sphere contribution l s arising from the change in the environment's polarizability and medium relaxation due to charge transfer: The intra-molecular contributions associated with the transfer of one electron were evaluated by using the adiabatic potential (AP) approach, which is also known as the "four point method".The outer sphere contribution cannot be so easily estimated for charge transport processes in solid state media.In the past few years, some efforts have been made in this direction by the groups of Brédas, 42 Troisi 43 and Andrienko, 44 who evaluated l s by hybrid QM/MM methods.However, only the class of acene compounds has been studied in depth.The use of continuum theory is another possible route, which has been followed by Gajdos et al. 28 According to their evaluations, l s in PCBM crystals turns out to be in the range of 25-36 meV.In our study l s has been taken equal to 36 meV, the upper limit of this range, thus enhancing charge localization.
Another important parameter in eqn ( 1) is the effective nuclear frequency.To access this parameter, the Huang-Rhys (HR) factors S j (electron-phonon coupling terms) have been computed from the displacement parameters (see the ESI †) as in ref. 45 and 46.Since low-frequency vibrations can be described to a good approximation in classical terms, and because of their possible anharmonicity 47 (they are typically associated with large-amplitude librations around single bonds), the contributions for vibrational modes below 250 cm À1 were not included in the evaluation of the effective nuclear frequencies of PCBM and ThCBM.This cutoff falls below the lowest vibrational frequency of C 60 , so that all vibrations associated with the fullerene cage are treated quantum mechanically, and all the low-frequency "classical" vibrations are associated with the side groups.The effective nuclear frequency can then be written as: where n j are the normal mode frequencies above the cutoff and S j are the associated HR factors.For the evaluation of l i and n eff , the hybrid B3LYP DFT functional and the double split 6-31G** basis set have been considered.The equilibrium geometries of the neutral and anionic states of the different species have been fully optimized, and the respective vibrational force elds have been evaluated to compute the HR factors.These calculations involved a standard, single-determinant representation of the electron density.In principle, a multi-conguration approach would have been more accurate (with an appropriate choice of the active space!) and consistent with our ET simulations including the low-lying L1 and L2 orbitals, but presently such calculations are still too demanding for molecules of this size.All these calculations have been carried out with the Gaussian09 program. 48

Results and discussion
Overall, we have investigated six crystal structures.Four are PCBM co-crystals obtained by crystallization from dichloromethane (DCM), 25 monochlorobenzene (MCB), 20 ortho-dichlorobenzene (DCB) 20 and CS 2 , 24 all of which contain guest solvent molecules.One is the solvent free structure (n-PCBM) recently described by our group 22 (related to the one obtained from DCB) and the last one is a co-crystal of ThCBM, recently obtained from CS 2 by Choi et al. 24 Table 1 summarizes the space group, the number of solvent molecule per unit cell as well as the PCBM (ThCBM) coordination numbers for each crystal structure (the latter can be dened as the number of fullerene molecules directly adjacent to a given one).Histograms (not shown) of the fullerene-fullerene distances were computed to determine a reliable cutoff (1.1 nm) for the denition of the rst nearest-neighbours for each crystal structure.Four of the structures are monoclinic (P2 1 /c or P2 1 /n space groups).Their unit cells contain four PCBM molecules, all of which are related by symmetry.The two remaining structures are triclinic (P 1 space group) and they contain two non-equivalent pairs of symmetry-related fullerene molecules.In the triclinic structures we report two coordination numbers, respectively for the rst and second sets of non-equivalent molecules, which in one case (ThCBM/CS 2 ) are different from each other.
In the following section we: (a) report the calculation of the reorganization energies and effective mode frequencies of PCBM and ThCBM and compare then with those of C 60 , (b) discuss the electronic couplings and site energies, including the effects of theory level and solvent molecules, and identify the conditions for the possible breakdown of the hopping model, (c) report and rationalize the mobility tensors corresponding to percolation networks of various complexities for each crystal structure.

Reorganization energies
The AP intra-molecular reorganization energy for C 60 (used as a reference system for comparison), PCBM and ThCBM are reported in Table 2. Values for C 60 and PCBM agree with literature data 27,30,49 while for ThCBM this is, to the best of the authors' knowledge, the rst evaluation of l i .C 60 features the lowest reorganization energy among all the compounds, since the extra electron is fully delocalized over the p-electron carbon cage.PCBM and ThCBM show slightly higher reorganization energies (>0.015 eV) due to the phenyl and thienyl lateral groups which change locally the p-electronic structure of the fullerene and introduce extra vibrational modes in the charge relaxation process.
To better understand the role of the local electron-phonon coupling parameters upon charge transfer, l i has been decomposed over all the vibrational normal modes of the neutral and charged species and both contributions (i.e.projection of the neutral geometry on the anion state and vice versa) are reported in Fig. 2 for each compound (numerical values are in the ESI †).We can observe that the introduction of the phenyl (PCBM) and thienyl (ThCBM) lateral group on the C 60 carbon cage perturbs the local electron-phonon couplings changing those normal modes sensitive to the electron transfer process.Overall, the values of S j terms are higher for PCBM or ThCBM with respect to the C 60 case.PCBM and ThCBM present a few normal modes very sensitive (i.e.high reorganization) to the charge transfer process, showing HR factors of the order of 0.1.In particular, for PCBM we have the 903 cm À1 (S ¼ 0.12), 1490 cm À1 (S ¼ 0.09) and 1800 cm À1 (S ¼ 0.06) active vibrations, while 902 cm À1 (S ¼ 0.28), 1490 cm À1 (S ¼ 0.095) and 1827 cm À1 (S ¼ 0.16) for ThCBM.These vibrations represent C]C/C-C oscillations of the C 60 cage and C]O stretching of the carbonyl bond.
On the basis of these HR factors we computed the n eff values compiled in Table 2.The high effective frequencies obtained for PCBM and ThCBM reect not only the higher reorganization energy with respect to C 60 but also the high anharmonicity of low energy normal modes introduced by the side groups.The l i and n eff values reported in this section have been used in eqn (1), with the other ET parameters, for the evaluation of the rates.

Electronic properties and adiabaticity
Electronic couplings were calculated at the DFT level using the B3LYP functional and the 3-21G basis set.This choice was made  This journal is © The Royal Society of Chemistry 2014 on the basis of preliminary calculations on seven different PCBM dimers.The couplings between their L0, L1 and L2 orbitals were computed at the DFT level using the B3LYP and BLYP functionals and the 6-311G**, 6-31G*, 3-21G and STO-3G basis sets, and at the semi-empirical level using the ZINDO/S, AM1 and PM3 parameterizations.Fig. 3 shows a correlation diagram of these electronic couplings which demonstrates the substantial agreement among the DFT calculations, with the exception of those using a minimal basis set.The latter provide couplings which are intermediate between the non-minimal DFT and the rescaled AM1 or PM3 calculations (the correlations represented in Fig. 3 are independent of such scaling).It is also interesting to point out that, even if the ZINDO/S electronic couplings are comparable in magnitude to the DFT calculations (this is well known, and is one reason why they have been frequently used 17,26,27 ), Fig. 3 indicates that ZINDO/S results do not correlate so well with DFT results.The scaling factor between the semi-empirical results and the DFT results obtained with the largest basis set (6-311G**) has been computed by linearly tting the data obtained at each level of calculation plotted as a function of each other (not shown).The scaling factors (the slope of the linear ts) reported in Table 3 conrm that the ZINDO/S electronic couplings are closer to the DFT than AM1 or PM3.
Site energies were estimated at the semiempirical level using the AM1 parameterization.This choice was dictated by the need to maintain reasonable computational times even for large clusters of fullerene and solvent molecules, which are responsible for sizeable polarization effects.A correlation diagram similar to Fig. 3 has been included in the ESI † (see Fig. S4) for the site energy differences.It demonstrates that, even though there are some differences among the computational levels, these are much more limited than for the couplings, since the correlation coefficients obtained from the dimer calculations are always comprised between 0.94 and 1.The scaling factors in Table 3 support our choice, as the AM1 parameterization is the one which provide the best tting of the DFT results.surrounded by their nearest neighbors, excluding the solvent molecules and (c) fullerenes with their nearest neighbors, including the solvent.In the latter case, when the solvent molecules have some positional/orientational disorder (PCBM/ DCM and PCBM/MCB), site energies were calculated for all possible congurations.Finally, for the sake of comparison, we also give (d) the site energies extracted from the dimer calculations (i.e., those used for the evaluation of the electronic couplings).
The position of the L0, L1 and L2 levels calculated on isolated molecules depends on the crystal structure because of subtle differences in the PCBM (ThCBM) molecular geometries.In particular, the conformation of the PCBM and ThCBM side chains can lead to specic intra-molecular interactions which shi these levels.For the same reason, non-equivalent molecules in triclinic structures also present slightly different L0, L1 and L2 energies.The L0L0, L1L1 and L2L2 site energy differences between non-equivalent sites may be neglected when considering isolated molecules, but they become much larger when including polarization effects from neighbouring fullerene derivatives (b).Thus, the L0 and L1 as well as the L1 and L2 levels of non-equivalent sites may become very close in energy and they actually cross each other in the case of the ThCBM/CS 2 structure (i.e., the L1 level of one molecule is lower than the L0 level of another, and similarly for the L2/L1 levels).The smallness of L0L1 and L1L2 site energy differences points to the possible contribution of the higher energy virtual orbitals to ET in the triclinic structures.At the same time, the site energies for the monoclinic structures suggest that only L0 contributes to ET.These indications have been checked by explicit KMC simulations, to be discussed below.
The inclusion of the solvent molecules in the coordination spheres (c) produces an additional shi of the fullerene levels.This shi appears to be more pronounced for the polar solvents.Indeed, when comparing the PCBM/DCM and PCBM/CS 2 structures, which belong to the same space group and have very similar lattice parameters, we nd that the levels of the PCBM/ CS 2 structure are almost unchanged while those of the PCBM/ DCM structure are slightly shied toward higher or lower energy, depending on the positions/orientations of the solvent molecules.This effect is even more pronounced for the PCBM/ MCB structure, as the levels obtained for the different positions/ orientations of the solvent molecules are spread to the point that they form a near-continuum structure.This is due partly to the polarity of MCB, partly to the fact that there are three solvent molecules surrounding each PCBM in this structure (DCM is also polar, but there are only two solvent molecules adjacent to the PCBM in the PCBM/DCM structure).
The calculations in case (c) show that, even though solvent molecules do not directly participate in charge transport (they have high energy orbitals and large reorganization energies), they might have an indirect effect by producing a spread in the fullerene orbital energies.However, the importance of this energy spread depends also on the solvent dynamics.For slow solvent reorientational motions, the electron would "see" a static distribution of energy levels whereas for fast solvent motions, it would feel an effective, narrower distribution centered about the average value.It turns out that the structure with the wider distribution (PCBM/MCB) is also the structure which, according to our previous MD simulations, 19 has the fastest solvent rotations.For this reason, and for the sake of simplicity, we have neglected the solvent contribution to the site energies and used those from the (b) case in our KMC simulations.
The results of (d) case calculations on the fullerene dimers are rather far both from the isolated molecules and from the molecules within a large cluster.This conrms that the site energies from the dimer calculations are ill-dened and they cannot be used reliably in a KMC simulation.
Fig. 5 contains level plots of DE ‡ ad as a function of the orbital couplings and energy differences, for each crystal structure (in principle, the barrier depends also on the reorganization energies, but these are xed parameters here).For some specic transport directions, we obtain negative adiabatic barriers for ET (points within the grey regions in the gure).This was observed also by Blumberger et al., 28 who justied this result on the basis of the inter-orbital couplings (for large couplings, an electron cannot be considered localized, as it spreads over strongly coupled molecules).However, these authors restricted their analysis to the L0L0 transfers and neglected the effect of site energy differences, even for the triclinic crystal structures with non-equivalent fullerene molecules.Our results for the sign of DE ‡ ad are somehow different from theirs, due to smaller electronic couplings and to the fact that our calculations include the effect of site energy differences.These cannot be neglected, both because we consider three orbitals instead of one and because some of our crystal structures (the triclinic ones) include PCBM or ThCBM molecules which are unrelated by symmetry.
Three different types of ET can be distinguished in Fig. 5: "at" (L0L0, L1L1, L2L2; green symbols), "upward" (L0L1, L1L2, L0L2; blue symbols) and "downward" (L1L0, L2L1, L2L0; red symbols).For the monoclinic structures, as no positional disorder has been introduced in our model, the at, upward and downward ETs are respective athermal (DE if ¼ 0), endothermal (DE if > 0) and exothermal (DE if < 0).For the triclinic structures, for which the L0, L1 and L2 levels of non-equivalent sites can cross each other, the sign of the site energy difference cannot be related any more to the type of ET (for example, there can be upward but exothermal L0L1 electron transfers).
It appears that for the monoclinic structures the ET pathways within the F1-networks dened in Section 2.2 (light green squares in the rst four panels) never reach the negative region of the maps.This suggests that the hopping model is acceptable for these percolation networks including only LUMO-LUMO transitions.However, looking at the F3-networks (all red, green and blue points), a major part of the exothermal pathways (red points) display negative DE ‡ ad .For the triclinic structures, even the F1-networks display some negative barriers (light green squares in the last two panels), corresponding to ET between non-equivalent sites with different LUMO energies.
To summarize, in our model DE ‡ ad can be negative for two distinct reasons: (i) large electronic couplings, leading to charge delocalization as proposed by Blumberger and coworkers, 28 but also (ii) large site energy differences, preventing charge localization on high energy sites even for fairly small electronic couplings.In principle, low energy fullerene sites may act as shallow traps driving ET.The KMC simulations will allow us to check this hypothesis.

Percolation networks
Using the electronic couplings, site energies and reorganization energies described in the previews sections, ET rates have been evaluated for each crystal structure according to the LZ expression (see the ESI † for the Marcus and MLJ rates).Rather than presenting large tables of numbers, which would be of limited interest due to the approximations in our calculations, we will resort to graphical representations produced with the "Rpdb" 51 and "rgl" 52 packages.These serve better our main purpose, that is, to get a feeling for the consequences of molecular packing on the ET properties of these fullerene derivatives.
Fig. 6 presents the ET rates which dene the F1-network for each crystal structure.The rates are represented by orange cones, whose tips point towards one molecule and whose bases are proportional to the rate to that molecule (i.e., the radius is proportional to the square root of the rate).The reader may think of them as funnels, channeling electrons towards their destination site.Note that the scaling factor connecting the size of the cones to the ET rates is the same of all crystal structures.Therefore, they provide a clear representation of the dimensionality of the percolation networks (one-, two-or threedimensional transport), of the main electron transport directions, and of the relative ET rates within different crystal structures.
The cones in Fig. 6 show that the ET pathways within PCBM/ DCM and PCBM/CS 2 form a well-developed three-dimensional network, with a certain preference for the a-axis.Note the strong similarity between transfer networks in these two crystals, which reects the analogy in their fullerene packings.Upon closer scrutiny, it turns out that the ET rates along a are actually 1.3 times larger in PCBM/DCM, in agreement with the slightly lower value of this lattice parameter.The two other monoclinic structures have lower dimensionality percolation networks.PCBM/DCB is essentially two-dimensional, in agreement with the alternation of PCBM and solvent layers along the (1 0 1) direction.Surprisingly, the n-PCBM network is the one with the lowest connectivity.The graphical representation of the rates indicates that the electrons may zig-zag comfortably only along the b-axis.A very short contact distance (2.9 Å) between PCBM molecules along that direction is at the origin of these zig-zags.Finally, the ET networks of the triclinic structures are more complex due to the presence of two distinct sets of non-equivalent fullerenes.The hopping rate from one to the other set can be one order of magnitude (or more) larger than the opposite one, and this is borne out by the very asymmetric double cones connecting them.It is difficult to decide on the dimensionality of their percolation networks by looking at these gures.In principle, the fullerenes at the end of the thickest cones might dominate charge transport.
The ET rates represented in Fig. 6 correspond to the L0L0 transitions forming the F1-networks.It is extremely difficult to obtain and interpret a similar representation of the F3networks, which include transitions to/from the L1 and L2 orbitals (there are actually 18 such transitions between any pair of neighbouring fullerenes!).To understand the role of the higher energy orbitals and to answer some of the questions raised above, we carried out the KMC simulations described in the following section.

Electron transport simulations
To conrm and extend our previous conclusions on the charge transport properties of the crystals, which were based on the rates for the L0L0 transfers, we carried out KMC simulations on the F1-and F3-networks.For the triclinic structures, the role of high energy molecules was also investigated by constrained KMC simulations on the R1-and R3-networks.The numerical results for the electron mobilities are collected in Table 4 and, for the F1-and F3-networks, they are also represented graphically in Fig. 7.
Fig. 7 illustrates both the mobility tensors (ellipsoids) and the recorded probabilities of the electron hops between neighbouring molecules in the F1-and F3-networks (respectively on the le and right part of the gure).For the probabilities, we use again a representation in terms of cones, but now the area at the base of each cone is proportional to the number of hops towards a given site.Note that the pair of cones connecting two fullerenes is symmetric, even in the triclinic structures with high-and low-energy molecules.The reason is the compensation between site-occupancies (P i ) and hopping probabilities (p ij ), which follows from the detailed balance condition: P i p ij ¼ P j p ji .
Table 4 Eigenvalues (cm 2 V À1 s À1 ) and eigenvectors of the mobility tensors calculated considering either the F1-, F3-, R1-or R3-networks (using the LZ rates) for each crystal structure We rst discuss the KMC results for the F1-networks, which neglect the effect of higher energy orbitals (Fig. 7 le).First of all, the dimensionality of the electron transport networks is conrmed for all monoclinic structures (3D for PCBM/DCM and PCBM/CS 2 , 2D for PCBM/DCB, 1D for n-PCBM).The simulations conrm also the similarity between the rst two structures.Their mobility ellipsoids are strongly oriented along the short a-axis, but they have non-negligible components also in the orthogonal directions.Like the hopping rates, the ratio of their long axes is close to 1.3 (see Table 4).The two other monoclinic structures forbid electron transport along one or two directions, as indicated by the presence of near-zero mobility components.Going to the triclinic structures, PCBM/ MCB is the one with the lowest mobilities.Transport is threedimensional, but even the largest component is relatively small.Strikingly, in this structure, both high-low-energy molecules contribute to transport, as conrmed by comparison with the reduced R1-network (the largest mobility is reduced by a factor of four, and its two orthogonal components become exactly zero).Instead, the F1 and R1 mobilities are very similar in the ThCBM/CS 2 structure.Thus, only low energy fullerene molecules contribute to charge transport, which is strongly one dimensional and oriented along the a-axis, in agreement with experiment. 24Here the "electron channels" are perfectly linear, unlike those of n-PCBM which appear to be zigzag-like.The large mobility associated with these channels (1.48 cm 2 V À1 s À1 ) is in good agreement with the experimental value (2 cm 2 V À1 s À1 ), considering the approximations in our electron hopping model and the neglect of thermal motions.We also found for the ThCBM crystal structure a very large anisotropic factor as experimentally observed by Fukuzumi et al. 24 Fig. 7 Effective hopping probabilities (cones) and associated mobility tensors (ellipsoids) within the F1-networks (left) and F3-networks (right), obtained for the KMC simulations with the LZ rates.
The analysis of the site energies reported above suggested that triclinic should be more sensitive than monoclinic structures to the inclusion of higher energy orbitals (F3-networks, see the right-hand-side of Fig. 7).Indeed, it has only little effect on the rst two monoclinic structures but it does have an effect on the other two.While the two non-zero mobility components of the PCBM/DCB structure are increased by 20-30% and transport remains two-dimensional, the n-PCBM structure switches from almost perfectly one-dimensional to a highly anisotropic three-dimensional transport.Similarly, the higher energy orbitals increase the average value as well as the isotropy of the mobility within the PCBM/MCB crystal.Instead, relatively minor effects are observed in the ThCBM structure which remains the most efficient but one-dimensional electron transporter.This is also conrmed by the R3-network simulations.
The individual ET rates can vary between 10 4 and 10 13 s À1 , with signicant differences among the different crystal structures.Their values are given in the form of histograms in Fig. 8.They are decomposed based on whether they involve only the F1-network (L0L0 transitions), they lead from the F1-network to some higher energy orbital (L0L1 and L0L2 transitions), or otherwise (L1L1, L2L2, L1L2, L2L1, L2L0, and L1L0 transitions; the ESI † contains further histograms with a complete decomposition by the type of transition).The structures which are insensitive to the L1 and L2 orbitals (PCBM/DCM and PCBM/ CS 2 ) have an order-of-magnitude separation between the L0L0 transitions and the L0L1 and L0L2 ones.Thus, the higher energy orbitals are irrelevant and transport is conned to the F1-network.In the other structures, there is some overlap between the distributions of these types of transitions, implying a competition between ground-and excited-state transport.The low sensitivity of the ThCBM/CS 2 structure despite this overlap is due to a particular connectivity of its networks which conne the transport to the R1-network.

Summary and conclusions
Progress in organic electronics and photovoltaics requires a thorough understanding of structure-property relationships.These include not just the molecular structure and properties, which from a computational point of view can be now tackled by standard quantum chemical methods, and also supramolecular aspects arising from intermolecular interactions.These are much harder to approach by purely computational methods and, as result, it is almost imperative to combine computational and experimental data on well-dened systems.We have taken such an approach in this study, which has focussed on the relationship between molecular structure, intermolecular packing and electron transport properties for PCBM-currently the most widely used n-type material-and ThCBM-a new, closely related fullerene derivative.
We have considered all the available experimental crystal structures of these compounds and applied a hopping model based on a Landau-Zener estimate of the charge transfer rates.Mobility tensors were evaluated by KMC simulations which, in some cases, included the possibility of electron transfer events to/from higher-than-LUMO orbitals.Our results can be summarized as follows: PCBM-like molecules may crystallize in different polymorphs forming one-, two-or three-dimensional percolation networks.
The dimensionality of the electron percolation networks depends on the distribution of the inter-fullerene distances.High electron mobility directions are typically associated with short inter-fullerene distances.For this reason, crystal structures with a short lattice vector display high electron transporting linear channels along that direction.
The fullerenes' LUMO + 1 and LUMO + 2 orbitals may participate in electron transport.This is especially important for triclinic structures with non-equivalent fullerene molecules as, for example, the LUMO + 1 orbital on one molecule may be degenerate with the LUMO on another.By extrapolation, we expect this "excited state conduction" to be even more important for amorphous structures, where each molecule has a different geometry and environment.
The highest mobility value has been obtained for a particular direction of the ThCBM/CS 2 structure, in which This journal is © The Royal Society of Chemistry 2014 electron transport is highly anisotropic because of perfectly aligned channels presenting high electron transfer rates.We nd a mobility of 1.48 cm 2 V À1 s À1 associated with that direction which is in good agreement with the experimental mobility value recently reported by Fukuzumi et al. (2.0 cm 2 V À1 s À1 ).ThCBM molecules which dominate charge transport are actually in a low coordination environment (coordination number ¼ 7), while there is a second set of molecules with a higher coordination (10) which make a minor contribution to the mobility.
Further work is currently in progress in order to identify and characterize other crystalline forms of PCBM and new fullerene derivatives. 25The availability of further structural and mobility data on well-dened crystal forms will provide a useful testing ground for current charge transport models.Indeed, there is a real need for new models and computational approaches, overcoming the picture of discrete hopping events among localized electronic states. 53,54Its fundamental inadequacy, which for PCBM already been pointed out by others, 27,28 emerges here from the fact that some of our calculated ET rates are comparable with the vibrational frequencies of the fullerene cage (10 13 s À1 x 300 cm À1 ).As a result, the time scale separation assumed by transition-state-type approaches simply does not apply.Suitable extensions of a recently proposed coarsegrained quantum-chemical model of these materials 55 may prove to be useful in this context.

Fig. 4
collects the site energies for all crystal structures, in the cases with (a) isolated fullerene molecules, (b) fullerenes

Fig. 3
Fig.3Correlation diagram (obtained with the "corrplot" library 50 of the R project 32 ) of the electronic couplings calculated between the L0, L1 and L2 orbitals of seven different PCBM dimers at different DFT and semi-empirical levels.

Fig. 4
Fig. 4 Site energies (L0 in red, L1 in green and L2 in blue) within the different crystal structures as calculated on (a) isolated fullerenes, (b) fullerenes surrounded by their first coordination spheres, removing solvent molecules, (c) fullerenes surrounded by their first coordination sphere, including solvent molecules and (d) fullerene dimers.

Fig. 5
Fig. 5 Level plot (in the [H if , DE if ] plane) of the adiabatic free energy barriers (in eV) for the different crystal structures.The grey area indicates the negative region of the plot.The green, red and blue symbols indicate respectively the flat, downward and upward ET directions.

Fig. 6
Fig.6Representation of the Landau-Zener ET rates for the L0L0 transitions (cones with bases proportional to the rates, pointing towards the destination fullerene), forming the F1-networks.

Fig. 8
Fig. 8 Decomposition of the LZ ET rate distributions of the different crystal structures.Green and blue bins are respectively the contributions of transitions forming the F1-network (L0L0) and those allowing escape from it (L0L1 and L0L2).

Table 1
Structural information of the different crystal structures

Table 2
AP intra-molecular reorganization energy and associated effective mode energy/frequency for C 60 , PCBM and ThCBM

Table 3
Scaling factors for the semi-empirical (AM1, PM3 and ZINDO/ S) results to produce electronic couplings (top) and site energy differences (bottom) fitting with the DFT (B3LYP and BLYP) results obtained with the largest basis set (6-311G**)