Linear-scaling density functional simulations of the eﬀect of crystallographic structure on the electronic and optical properties of fullerene solvates †

In this work, the crystal properties, HOMO and LUMO energies, band gaps, density of states, as well as the optical absorption spectra of fullerene C 60 and its derivative phenyl-C 61 -butyric-acid-methyl-ester (PCBM) co-crystallised with various solvents such as benzene, biphenyl, cyclohexane, and chloro-benzene were investigated computationally using linear-scaling density functional theory with plane waves as implemented in the ONETEP program. Such solvates are useful materials as electron acceptors for organic photovoltaic (OPV) devices. We found that the fullerene parts contained in the solvates are unstable without solvents, and the interactions between fullerene and solvent molecules in C 60 and PCBM solvates make a significant contribution to the cohesive energies of solvates, indicating that solvent molecules are essential to keep C 60 and PCBM solvates stable. Both the band gap ( E g ) and the HOMO and LUMO states of C 60 and PCBM solvates are mainly determined by the fullerene parts contained in solvates. Chlorobenzene-and ortho -dichlorobenzene-solvated PCBM are the most promising electron-accepting materials among these solvates for increasing the driving force for charge separation in OPVs due to their relatively high LUMO energies. The UV-Vis absorption spectra of solvent-free C 60 and PCBM crystals in the present work are similar to those of C 60 and PCBM thin films shown in the literature. Changes in the absorption spectra of C 60 solvates relative to the solvent-free C 60 crystal are more significant than those of PCBM solvates due to the weaker eﬀect of solvents on the p -stacking interactions between fullerene molecules in the latter solvates. The main absorptions for all C 60 and PCBM crystals are located in the ultraviolet (UV) region.


Introduction
Conjugated polymer/fullerene organic solar cells (OSCs) have high potential as low-cost but efficient renewable energy sources and have attracted increasing attention in recent years, as they are relatively lighter, cheaper and easier to produce from renewable sources in contrast with most inorganic solar cells. [1][2][3][4][5][6] Today, power conversion efficiencies (PCEs) of over 10% for solar cells based on blend films of conjugated polymers with fullerenes have been reached. 7,8 In these high performance devices, fullerenes are used as electron-accepting materials. Furthermore, it has been found that fullerenes crystallize in some solid blends prepared from organic solvents such as benzene, chlorobenzene (CB), ortho-dichlorobenzene (oDCB), and chloroform at high concentration or following annealing treatment. [9][10][11][12][13][14][15][16][17][18][19][20][21][22][23] The crystallite formation of fullerenes can have a significant effect on the properties of blend films because the properties of crystalline fullerenes may be different from those of amorphous fullerene molecules. For instance, the band gaps (E g ) of both C 60 and its derivative phenyl-C 61 -butyric-acid-methyl-ester (PCBM) crystals are lower when compared to those of the isolated C 60 and PCBM molecules. 24 This E g reduction has a significant effect on their optical properties. Moreover, the large crystallite formation of fullerenes in some systems may aid charge mobility relative to the amorphous phase due to the molecular ordering and the lack of grain boundaries. 25,26 Some previous studies 14,16,17,20,21,23 have also shown that devices based on poly-3-hexylthiophene (P3HT) and PCBM blends display a better performance (e.g., a large increase in photocurrent and carrier lifetime) following thermal annealing which alters the fraction of amorphous and microcrystalline PCBM domains, with the formation of PCBM nano-crystalline domains that are crucial for high PCEs. 14,16,27 To date, a large number of studies both in experimental implementations [9][10][11][12][13][14][15][16][17][18][19][20][21][22][23][27][28][29][30][31][32][33][34][35][36][37][38][39][40][41][42][43] and in theoretical aspects 24,[44][45][46][47][48][49] have been carried out on the preparation, characterization of morphology and structure and physical properties of solvent-free or solvated C 60 and PCBM crystals. All these studies can not only help to improve the performance of these crystals in OSCs but can also help to provide a deeper insight into their fundamental properties. It is found that both solvent-free C 60 and PCBM crystals crystallize with four molecules, in space group Pa% 3 of the cubic system 34,35 and P2 1 /c of the monoclinic system 37,38 respectively [see Fig. 1(a) and (e)], while the crystal structure of solvated C 60 and PCBM co-crystals strongly depends on the solvent because it can diffuse into the crystal lattice and play an important role in modifying the interactions between the molecules; thus fullerenes can assume different arrangements when using different solvents. For instance, Meidine et al. 33 found that the crystal structure of benzene-solvated C 60 is triclinic with a P% 1 space group symmetry where the ordered molecules of C 60 are in an approximately hexagonal close-packed arrangement and separated by benzene molecules. Pénicaud et al. 30 yielded some single crystals of the formula [(C 60 )][(C 6 H 5 ) 2 ] with a monoclinic structure by slow co-crystallization of saturated solutions of C 60 and biphenyl in toluene at ambient temperature. Rispens et al. 27 obtained two PCBM packing structures from single crystals grown from CB and oDCB by using extremely slow solvent evaporation conditions and maintaining a solvent saturated atmosphere. When PCBM crystallizes from CB, the co-crystal structure is triclinic with four PCBM and two solvent molecules in each unit cell, while a monoclinic structure with four PCBM and four solvent molecules in the unit cell is formed when using oDCB (the asymmetric units of these crystal structures are shown in Fig. 1). From a structure-property relationship point of view, different crystal structures may result in different electronic and optical properties. Both of these are very important for the optimization of solar cell devices, in particular since the maximum open-circuit voltage (V oc ), one of the most important parameters for the solar cell efficiency, largely depends on the gap between the highest occupied molecular orbital (HOMO) energy of the donor material and the lowest unoccupied molecular orbital (LUMO) of the acceptor material. 50,51 In addition, the photocurrent is also determined to some extent by the absorption coefficient of the photoactive layer. 52 In contrast to the case of solvent-free C 60 and PCBM crystals, whose electronic and optical properties have been widely investigated, 31,[40][41][42][43][44][45][46][47][48][49] there is little information available on the electronic and optical properties of solvated fullerene crystals so far, to the best of our knowledge, and therefore, in order to tailor the efficiency of polymer/fullerene based solar cells in a suitable way for higher PCEs, it is necessary to study the electronic and optical properties of fullerenes co-crystallised with different solvents. To this end, a suitable approach is to use simulations with first-principles density functional theory (DFT) that can provide accurate ground state properties such as energy eigenstates and also allows the calculation of optical absorption spectra. Conventional first-principles DFT codes such as CASTEP, VASP, PWSCF and ABINIT require a computational effort that scales cubically with the model size. As the solvated co-crystals we are interested in here are large in scale, involving as many as 216 to 2208 atoms, calculations using these codes can be too computationally expensive. Hence, we have employed a linear-scaling (LS) DFT method as implemented in the ONETEP (Order-N Electronic Total Energy Package) program 53 which is particularly suitable for studying large systems, since the computational cost only increases linearly with respect to the number of atoms.

Computational details
All structural properties of fullerene crystals such as the space group and lattice parameters in this paper (summarised in Table 1) are taken from the Cambridge Structural Database (CSD), except for the 1-methylnaphthalene-solvated PCBM co-crystal whose crystallographic structure is reported herein for the first time (see ESI †). One should note that these structures may not be the only possible ones since fullerene crystals usually can have a number of polymorphs close in energy, and there is no guarantee that these polymorphs whose X-ray structures are solved in bulk powder samples are prevalent or even present in films. Although this is not the focus of this work, in order to predict the possible or the most stable structures of these fullerene crystals, computational methods for organic crystal structure prediction based on global minimisation of the lattice energy 54 can be employed, since they have been successfully applied to the prediction of single-55 and multicomponent 56 organic molecular crystals and their reliability has been tremendously improved in recent years. 57 However, XRD is the most reliable method for the experimental determination of crystal structures and the R 1 and wR 2 refinement indices for the structures we considered are within accepted limits.
During the geometry optimization we kept the dimensions of the unit cell as well as the atomic positions fixed, with the only exception of the hydrogen atoms, since it is not usually possible to accurately determine in experiments their position.
As is well known, the virtual orbitals in the exact Kohn-Sham (KS) model of DFT are genuinely bound one-electron states at exactly the same (local) potential as the occupied orbitals, and represent excited electrons in the non-interacting KS reference system. 58 Their energies, under certain conditions, can be interpreted as ionization energies from the corresponding states since the occupied KS orbital energies have already been interpreted as approximate ionization energies. 59 As a result, we can easily interpret the excitation spectrum using the one-electron KS model. In contrast, according to the canonical Koopmans' theorem, the energies of virtual orbitals in the Hartree-Fock (HF) model are interpreted as approximate electron affinities to an extra electron based on the frozen-orbital HF approximation. 60 Moreover, in many cases the HF virtual orbitals are not bound states and are shifted to much higher energies than the KS virtual orbitals, with their shapes being very diffuse. 58,60 Therefore, the HF virtual orbitals are often unphysical. The ONETEP program, which employs the one-particle density matrix reformulation of KS DFT, can describe the occupied (valence) KS states very well as compared to the conventional DFT implementations. However, it is unable to represent the unoccupied (conduction) equally well because the local orbitals used in ONETEP are specifically optimised to describe the valence states. In ONETEP calculations, these localized orbitals are represented by a set of non-orthogonal generalized Wannier functions (NGWFs) 61 which are atom centered and strictly localized within a set radius. The NGWFs are represented by a basis set of periodic cardinal sinc (psinc) functions, 62 which is equivalent to a plane-wave basis set, and Table 1 Refcode (CSD entry identifier), composition, space group, lattice parameters and density (r) of the C 60 and PCBM crystals as well as the temperature (T) at which these parameters are determined

Refcode
Composition (per unit cell) Space group is optimized self-consistently during the calculations. Due to the self-consistent optimisation of the NGWFs and the psinc basis set, the distinguishing feature of ONETEP is that it is able to perform linear-scaling DFT calculations 63 with near-complete basis set accuracy, as is possible in conventional cubic-scaling DFT approaches. Typically, the unoccupied states that lie 1-2 eV above the Fermi level obtained from the NGWFs in a standard ONETEP calculation are higher in energy than those obtained from a conventional cubic-scaling DFT code such as CASTEP, 64 and some unoccupied states are even lost. 65 Therefore, a recently developed methodology in which a second set of localized orbitals (conduction NGWFs) is optimized to accurately describe the unoccupied states has been employed. Further details of this new methodology have been given elsewhere. [65][66][67] Once the optimization of the set of conduction NGWFs has finished, the valence and conduction NGWF basis sets can be combined into a new joint basis set which is able to well represent both the occupied and unoccupied KS states of the system. The optical matrix elements of the electronic transitions between valence and conduction states can then be calculated using this new joint basis using Fermi's golden rule, which has been implemented in ONETEP. The electronic and optical properties obtained using the new joint basis set are expected to agree well with those obtained using the conventional cubic-scaling DFT formulation.
Using the dipole approximation in which the exponential is expanded in a Taylor series and only first order terms are included, the imaginary part e 2 (o) of the frequency-dependent dielectric function is defined as, where O is the cell volume, the indices v and c refer to valence and conduction bands respectively,|c n k i and E n k are the nth eigenstate and corresponding energy at a given k-point, and q is the polarization direction of the photon with an energy of ho. Using ONETEP the imaginary part of the dielectric function of many systems such as the metal-free phthalocyanine and C 60 -conjugated polymer hybrids has been successfully calculated. 65,66 One can then also obtain the real part e 1 (o) of the dielectric function from e 2 (o) using the Kramers-Kronig relation as follows, where P denotes the principal value, and it is generally set to 1 in calculations.
The optical absorption coefficient a(o) indicates how far light of a particular wavelength can penetrate into a material before it is absorbed. It can be calculated directly from the dielectric function by using the following formula: where c is the speed of light.
All calculations, including geometry optimizations, ground states and conduction states optimizations, have been performed using a dispersion-corrected PBE exchange-correlation (PBE-D2) functional which is based on the correction according to Grimme,68 with G point sampling only and periodic boundary conditions. The reason why Grimme's DFT-D2 was used is that D2 dispersion correction, which only includes two-body energies and neglects many-body terms and any effects originating from an atom's environment, has been widely used due to its very low computational cost. 69,70 In addition, since many results in the available literature have also been obtained by using PBE-D2, 47 the use of the same dispersion-corrected exchange-correlation functional would make our calculated results more comparable to them. Moreover, the D3 correction is not currently implemented in ONETEP.
The projector augmented wave (PAW) method was used in geometry optimizations due to its many advantages such as the higher computational efficiency, all-electron precision, and a lower psinc kinetic energy cutoff to converge compared with norm-conserving pseudopotentials (NCPPs). However, the NCPPs were employed in ground-state and conduction calculations because the PAW method is not currently available for computing e 2 (o) during conduction calculations in ONETEP. The total energies were converged with respect to the psinc cut-off energy and NGWF radii. The psinc grid spacing was set to be equivalent to a kinetic energy cutoff of 800 eV for geometry optimizations and 1000 eV for ground-state and conduction calculations and no truncation was applied to the density kernel. Both geometry optimizations and ground-state calculations were performed using 1 NGWF per H atom and 4 NGWFs each per C, O and Cl atom with a fixed radius of 8.0 Bohr. Conduction calculations were performed with 5 NGWFs per H atom and 13 NGWFs each per C, O and Cl atom, and a fixed radius of 9.0 Bohr. The number of conduction states explicitly optimized was increased linearly with the number of atoms in the system. A Gaussian smearing width of 0.1 eV was used for the calculations of density of states (DOS) and for the dielectric function. The imaginary part of the dielectric function was calculated using the momentum operator formalism since the position operator is ill-defined under periodic boundary conditions. Growth of single crystals of PCBM was attempted using slow crystallization conditions to maximize our chances of success. We employed slow vapour diffusion of isopropanol into a concentrated solution of 1-methylnaphthelene (1-MN). Various concentrations of PCBM in 1-methylnapthalene were explored in order to balance the rate of crystal formation to achieve diffractable single crystals. Attempts to grow single crystals from tetralin and ortho-xylene did not yield success.
Non-planar molecules are known to promote solvated structures. [71][72][73] This is presumably due to the significant void space commonly observed in crystals formed from non-planar molecules. The void volume and void fraction were calculated for the studied fullerene crystals using Crystal Explorer 3.1 (isovalue: 0.002, quality: high), as shown in Table 2. In each case, lattice inclusion of solvates reduced the void fraction of the crystals in comparison to the non-solvated crystals, the only exception being cyclohexane. However, a closer inspection of the cyclohexane solvate crystal revealed that the increased void space in this case is due to the inefficient molecular packing of the cyclohexane layer formed between fullerene layers within the crystal and not due to an increased void space within the fullerene layers (See Fig. S1, ESI †).
Void filling may not be the only driving force for solvate inclusion within the crystal lattice. It appears that solventsolute interactions are prevalent in the crystal structure of both C 60 and PCBM in combination with aromatic solvents. The host-guest inclusions of aromatic molecules display close contact between the p centroid of solvents and the fullerene cage (pÁ Á ÁC 60 ). Fullerene, being an electron acceptor and p-acidic, demonstrates a closer molecular contact to p-basic solvents (1-MN, benzene) than p-acidic solvents (oDCB). A deeper understanding of the strength of these host-guest interactions can be explained by studying their cohesive energies.
Cohesive energies of solvent-free and solvated C 60 and PCBM crystals. As shown in Table 3, the cohesive energy E coh (total) of a fullerene crystal was calculated as E coh (total) = À(1/m)[E cryst À mE mol (solvent) À nE mol (fullerene)], where E cryst is the total energy per unit cell of a fullerene crystal under periodic boundary conditions, E mol (solvent) and E mol (fullerene) are the energies of an isolated solvent and fullerene molecules, respectively, and m and n are the number of solvent and fullerene molecules contained in an unit cell. In order to study the contributions of the fullerene-fullerene, solvent-solvent and fullerene-solvent interactions to the cohesive energies of fullerene crystals, we have also calculated the cohesive energies of the fullerene and solvent parts in fullerene crystals as they represent the fullerene-fullerene and solvent-solvent interactions, respectively. Here, the E cryst of the fullerene (solvent) parts were calculated by removing the solvent (fullerene) molecules contained in the unit cells of crystals and by keeping the same lattice constants and atomic positions as their parent unit cells. The contributions of the fullerene-solvent interactions to the cohesive energies of fullerene crystals were calculated as E coh (fullerene-solvent) = E coh (total) À E coh (fullerene-fullerene) À E coh (solvent-solvent).
From Table 3 it is possible to observe that the cohesive energies E coh (fullerene-fullerene) of the solvent-free C 60 crystal and PCBM parts in CB-and oDCB-solvated PCBM co-crystals compare quite well with the experimental and theoretical values in the literature. 46,47,74 The E coh (fullerene-fullerene) of the PCBM part in the CB-solvated PCBM co-crystal is higher than the corresponding one in the oDCB-solvated PCBM co-crystal, which is also in good agreement with previous calculations using the PBE-D2 functional. 47 We find that the E coh (fullerene-fullerene) of pure fullerene parts increases with the decrease in the cell volume. As the solvent-free C 60 and PCBM crystals do not contain any solvent molecules, they are closest-packed and their cell volumes are the smallest among pure C 60 and PCBM parts, respectively. Their cohesive energies are the highest because the cell volume (per fullerene molecule) can represent to some extent the average distance between fullerene molecules under periodic boundary conditions for these pure fullerene parts, and therefore the average distance will decrease as the cell volume per fullerene molecule decreases, while the interactions between fullerene molecules become stronger. The cohesive energy, which is the energy required to separate a solid into its elementary ''building blocks'', can become higher within a certain distance range, and this can also explain the computational results reported by Gajdos et al.: 47  From the cohesive energies E coh (fullerene-fullerene) in Table 3 we can also conclude that the pure fullerene parts originated from fullerene solvates are more likely to decompose due to their lower cohesive energies, in contrast with the solvent-free fullerene crystals that are able to exist stably under realistic conditions; thus these pure fullerene parts would be metastable or unstable, while the fullerene solvates are stable as observed in experiments. This indicates that the solvents contained in the fullerene solvates are essential for the stability of these solvates. Due to the existence of solvent molecules, the p-stacking interactions between fullerene and solvent molecules can be formed in the fullerene solvates. From the cohesive energies E coh (total), E coh (fullerene-fullerene) and E coh (solvent-solvent) in Table 3 it can be seen that the cohesive energies E coh (total) of these C 60 and PCBM crystals are mainly determined by the fullerene-fullerene and fullerene-solvent interactions therein, except for the cyclohexane-solvated C 60 whose E coh (total) is mainly determined by the interactions between solvent-solvent molecules and fullerene-solvent molecules, as the number of solvent molecules is much higher than that in other fullerene solvates. To sum up, the contributions of the interactions between fullerene-solvent molecules to the cohesive energies E coh (total) of C 60 and PCBM solvates are significant in particular for the benzene-and biphenyl-solvated C 60 , and this implies that solvent molecules are indispensable to fullerene solvates.

Electronic properties
HOMO, LUMO energies and energy gaps of solvent-free and solvated fullerene crystals. From Fig. 2 and the densities of C 60 and PCBM crystals shown in Table 1, we found that both the HOMO and LUMO energies for C 60 and PCBM crystals shift upward as the crystal density increases, which is in agreement with the experimental results that suggest that higher HOMO and LUMO levels of the PCBM film could be linked with the increased film density. 41 As can be seen in Fig. 2, both the HOMO and LUMO energies of C 60 solvates are lower than those of the solvent-free C 60 crystal, and a similar behaviour can also be seen in the 1-methylnaphthalene-solvated PCBM co-crystal; however the HOMO and LUMO energies of CB-and oDCBsolvated PCBM co-crystals are higher compared to the solventfree PCBM crystal energies. As a result, CB-and oDCB-solvated PCBM co-crystals should be the most suitable electron acceptor materials to increase the driving force of charge separation in solar cells due to their higher LUMO energies. In addition, the magnitude of the changes in the HOMO and LUMO energies of C 60 solvates relative to the solvent-free C 60 crystal is clearly larger than that of PCBM solvates compared to the solvent-free PCBM crystal. This is probably because the p-stacking interactions between C 60 cages in C 60 solvates are largely reduced relative to the strong p-stacking interactions in the pristine C 60 crystal, as the existence of solvents therein separates the C 60 cages and enlarges the distances between the cages. As shown in Table 3, this kind of interaction only takes up a small portion (less than 50%) of all the interactions in C 60 solvates. Table 3 Cell volumes and cohesive energies E coh (total), E coh (fullerene-fullerene) and E coh (solvent-solvent) for C 60 and PCBM crystals, pure fullerene and solvent parts contained in fullerene crystals compared with the available experimental and theoretical results in parentheses; values in brackets are percentage ratios of E coh (fullerene-fullerene), E coh (solvent-solvent) and E coh (fullerene-solvent) to E coh (total)  In contrast, the solvents contained in PCBM solvates do not seem to significantly reduce the p-stacking interactions between PCBM molecules, with these being the strongest among the interactions in PCBM solvates (see Table 3). As shown in Fig. 3, the calculated band gap (E g ) of the pristine C 60 crystal is 1.36 eV, and despite being underestimated compared to the experimental value of 2.3 eV, 75 it is close to the published theoretical values of 1.09 eV obtained using PBE 44 and 1.5 eV obtained using LDA functionals 46 . In addition, the E g of 1.28 eV for monoclinic solvent-free PCBM crystal found in this work is also comparable with those of 1.21 and 1.32 eV for simple cubic and body-centered-cubic PCBM crystals, both calculated by using the LDA functional in ref. 24. Hybrid functionals such as B3LYP, 76,77 which includes a fraction of exact Hartree-Fock exchange, would provide an even more accurate prediction of band gaps, with values closer to the corresponding experimental ones. However, hybrid functionals are not currently implemented for periodic systems in ONETEP. 78 Both the pristine C 60 and PCBM crystals have a lower E g compared to the isolated C 60 and PCBM molecules, and this means that E g decreases for C 60 and PCBM from the gas phase to solid phase. The E g for solvated C 60 and PCBM co-crystals are different from those of the solvent-free C 60 and PCBM crystals.
To shed light on the effect of fullerenes and solvents contained in the fullerene solvates on the E g of solvates, we then calculated the E g of the solvent and fullerene parts of the solvates, as shown in Fig. 3. Both the pure solvent and fullerene parts have the same lattice constants as their parent solvated co-crystals, as they are obtained by deleting the fullerene and solvent parts of the co-crystals, respectively. From Fig. 3 it is possible to observe that the E g of pure fullerene parts derived from solvates are also lower than the E g of the corresponding isolated fullerene molecules, except for the one derived from benzene-solvated C 60 whose E g is equal to that of the isolated C 60 molecule. By comparing the E g of fullerene solvates with those of the corresponding pure solvent and fullerene parts, we found that the E g of fullerene solvates are slightly smaller than those of the corresponding pure fullerene parts, this probably being due to the effect of the solvents contained in the solvates.
The reductions of E g of the biphenyl-solvated C 60 and 1-methylnaphthalene-solvated PCBM co-crystals are the largest among the C 60 and PCBM solvates, respectively. However, the difference in the band gap between the fullerene solvates and their corresponding fullerene parts is so small that the E g of fullerene solvates is nearly the same as that of the corresponding fullerene parts. In contrast, the difference in E g between solvates and their corresponding pure solvent parts is significantly larger: this implies that the E g of solvated fullerene co-crystals is mainly determined by the fullerene part contained therein. Furthermore, the difference in E g among solvated fullerene co-crystals should be attributed to the different E g of fullerene parts. The reason why the fullerene parts in solvates have different E g is that their packing structures are different. In general, the packing structure of fullerene molecules in fullerene solvates is affected by the solvent: the solvent Fig. 3 Energy gaps of isolated C 60 and PCBM molecules, solvent-free C 60 and PCBM crystals as well as solvated fullerene co-crystals and their corresponding fullerene and solvent crystals which only retain the fullerene and the solvent parts of the co-crystal via removing the solvent and fullerene parts, respectively. molecules with which the fullerene molecules co-crystallise play an important role in modifying the interaction between fullerene cages via forming p-stacking interactions with fullerene molecules. Thus, the packing structure of fullerene molecules in fullerene solvates is a result of the interaction between solvent and fullerene molecules. Different solvents may result in different fullerene packing structures, which can consequently lead to variations in E g of fullerene solvates.
Electronic density of states of fullerene crystals. Fig. 4 shows the comparison of the occupied parts of the total density of states (TDOS) for the pristine PCBM crystal and isolated PCBM molecule with the experimental and theoretical data in ref. 40. As the TDOS of the PCBM film measured by photoelectron spectroscopy (PES) in ref. 40 gave the binding energies as positive numbers, we have changed the energy of the TDOS to be positive to match the results of PES, and also shifted the energy to align the lowest energy peak of TDOS. In Fig. 4, one should note that 0 eV on the energy scale is the Fermi level, not the vacuum level since photoemission spectra are usually displayed in binding energies relative to the Fermi edge of the sample. As can be seen, the TDOSs calculated in this work agree well with the experimental photoelectron spectrum and Fig. 4 Comparison of the occupied parts of the total density of states (TDOS) for the solvent-free PCBM crystal and isolated PCBM molecule between our results and the available experimental photoelectron spectrum and theoretical data in the literature. The TDOSs reported here were generated with a bigger Gaussian smearing width of 0.2 eV to make them more comparable with the results in the literature. theoretical results which were calculated by using the DMol 3 package using the BLYP functional, thus proving that our results are reliable.
In Fig. 5, we only show the DOS of fullerene crystals near the Fermi level because it is the most important part in the DOS and governs the properties of many materials. From the total and local DOS it is clear that the LUMO states of all fullerene solvates are entirely derived from the fullerenes C 60 and PCBM, but the HOMO states of these solvates are not. The contributions of biphenyl and 1-methylnaphthalene to the HOMO states of their fullerene solvates are obvious, thus also indicating that the effects of biphenyl and 1-methylnaphthalene on the HOMO states of their fullerene solvates are distinctly stronger than the others. Similar effects of these two solvents on the E g of their fullerene solvates have also been previously found, however the contributions of these two solvents are minor: the HOMO states of their fullerene solvates are still mainly determined by the fullerenes. Therefore, we can conclude that fullerenes C 60 and PCBM contribute greatly to the HOMO and LUMO states of their solvates, which is in agreement with the results we found above: the E g of fullerene solvates are mainly determined by the fullerene parts contained in solvates.

Absorption spectra of solvent-free and solvated fullerene crystals
The ultraviolet-visible (UV-Vis) absorption spectra of solvent-free and solvated C 60 and PCBM crystals are shown in Fig. 6. From Fig. 6(a) it is possible to see that the absorption peaks of the solvent-free C 60 crystal occurring at 449.7 nm, 394.0 nm and 334.5 nm shift about 104 nm or more towards longer wavelengths as compared to the absorption peaks of the C 60 thin film (346 nm, 268 nm and 220 nm) measured by Zhou et al.; 31 this red-shift is most probably caused by the GGA-PBE functional which usually underestimates the energy gap. As the intensity of the absorption spectrum of the C 60 thin film is dimensionless, we are not able to make a quantitative comparison with that of the absorption spectrum of the pristine C 60 crystal. In Fig. 6(e), a similar red-shift of the absorption peak (at 417.3 nm) of the solvent-free PCBM crystal can also be found due to the use of the GGA-PBE functional when compared to the absorption spectrum of the PCBM thin film, as recorded by Cook et al. 43 However, the general trend of the absorption spectra is not affected by the systematic red-shift error introduced by the PBE functional, and therefore for pristine C 60 and PCBM crystals, the spectra in the present work are similar to those for C 60 and PCBM thin films shown in ref. 31 and 43. As shown in Fig. 6(e), a weak tail Fig. 6 Ultraviolet-visible (UV-Vis) absorption spectra of pristine and solvated C 60 and PCBM crystals compared with experimental data 31,43 as well as the contributions of solvent molecules contained in fullerene solvates to the total absorption spectra. extending as far as 900 nm can be found in the absorption spectrum of both the pristine PCBM crystal and thin film. This absorption tail has also been observed in the absorption spectra of the pristine C 60 crystal and thin film at long wavelengths [see Fig. 6(a)], which confirms the results obtained by Cook et al.: the origin of this absorption tail for the PCBM thin film is from its parent C 60 cages rather than the light scattering or other experimental artifacts. 43 In Fig. 6(a) the weak absorption in the region from 500 to 900 nm for the pristine C 60 crystal and thin film is due to the high number of optical transitions that occur in this region for the C 60 molecule: for instance, nearly all transitions to the first excited state such as the HOMO-LUMO transition are symmetrically forbidden because of the high symmetry of the C 60 molecule. Only transitions to higher energy bands such as the second and third lowermost conduction bands that lie above the LUMO band, which corresponds to the absorptions of shorter waves, can be allowed.
As shown in Fig. 6, the main absorptions for all fullerene crystals studied here locate in the UV region. Comparing the absorption coefficients of C 60 solvates with those of the solventfree C 60 crystals, we found the absorption peaks of the former are lower than the corresponding ones of the latter, with the exception of the enhanced absorption peak of benzene-solvated C 60 in the middle UV (MUV) region (200 to 250 nm) and the one of biphenyl-solvated C 60 near 200 nm. The enhanced MUV absorptions for benzene-and biphenyl-solvated C 60 co-crystals should be attributed to the contributions of the solvent molecules contained in these two solvates, as shown in Fig. 6(b) and (c). From Fig. 6(d) it is possible to see that there is no contribution from the cyclohexane molecules to the absorption spectrum of the corresponding C 60 solvate. The absorption coefficients of benzene-, biphenyl-and cyclohexane-solvated C 60 co-crystals in the visible region (450 to 570 nm) are much smaller than those of the solvent-free C 60 crystal, thus indicating that these three solvated C 60 co-crystals have a rather weak ability to absorb light in this region. All of the four C 60 crystals show no absorption when the wavelength is above 600 nm. In Fig. 6(f)-(h), similar reductions can also be found for the absorption peaks of PCBM solvates in the region from 300 to 900 nm when compared to the solvent-free PCBM crystal. However, their reductions are small relative to the ones of C 60 solvates, since the decrease of p-stacking interactions between PCBM molecules caused by solvents in PCBM solvates is smaller than that in C 60 solvates. In particular, the CB-solvated PCBM co-crystal has a very similar absorption spectrum to the solvent-free PCBM crystal, which implies that CB is the optimal solvent to keep the absorption property of the PCBM crystal at its strongest. The absorption coefficients of the three PCBM solvates in the visible region (420 to 760 nm) are almost the same. The difference in computed absorption coefficients between CB-and oDCB-solvated PCBM co-crystals is so small that it would be difficult to observe experimentally. From Fig. 6(f)-(h), it can also be seen that the absorptions of PCBM solvates around 250 nm increase because of the significant contributions of the solvents in this region.

Conclusions
In conclusion, linear-scaling density functional theory (LS-DFT) calculations, which allowed us to simulate large systems with thousands of atoms, have been carried out in order to clarify the crystallographic, electronic and optical properties of the electron-acceptor materials for organic photovoltaics (OPVs) C 60 and its derivative PCBM crystals, as well as their solvates co-crystallised with different solvents, also including the PCBMcyclohexane solvate composed of 2208 atoms in its unit cell. Our results can be summarized as follows: (i) The relatively low cohesive energies of pure C 60 and PCBM parts which originate from the C 60 and PCBM solvates compared with those of the solvent-free C 60 and PCBM crystals with closest-packing structures indicate that the fullerene parts are unstable and prone to decompose without solvents. The significant contributions of the interactions between fullerene and solvent molecules to the cohesive energies of C 60 and PCBM solvates also imply that solvent molecules are essential for C60 and PCBM solvates.
(ii) With an increase in the crystal density, both the HOMO and LUMO energy levels for C 60 and PCBM crystals move upward, in agreement with experimental observations. We have examined a variety of solvents and found that chlorobenzene and ortho-dichlorobenzene are the best from our set of solvents for the preparation of OPVs due to the relatively high LUMO energies of their PCBM solvates, which are helpful to increase the driving force for charge separation in OPV devices.
(iii) Both the energy gap (E g ) and the HOMO and LUMO states of C 60 and PCBM solvates are mainly determined by the fullerene parts contained in solvates. The solvents play an important role in determining the packing structures of fullerenes contained in C 60 and PCBM solvates, which lead to the different E g of C 60 and PCBM solvates. Thus, the effect of solvents on the E g of fullerene solvates is not due to the p-stacking interactions between solvent and fullerene molecules but merely due to the different fullerene packing structures in solvates.
(iv) In contrast to the significant changes in the intensity of the UV-vis absorption spectra of C 60 solvates as compared with the one of the solvent-free C 60 crystal, the intensity of the absorption spectra of PCBM solvates changes slightly, relative to that of the solvent-free PCBM crystal. This is because the reduction of p-stacking interactions between fullerene molecules caused by solvents in PCBM solvates is not as strong as in the case of C 60 solvates. The main absorption peaks for these fullerene crystals are located in the UV region. The contributions of solvent molecules to the absorption spectra of their respective fullerene solvates mainly occur in the middle UV region (200 to 300 nm) with the exception of cyclohexane whose contribution cannot be found in the entire absorption spectrum of the C 60 solvate.
This work was motivated by the need to examine the properties of solvated C 60 and PCBM co-crystals, as it is possible that they can be used instead of solvent-free fullerene crystals. We expect our obtained results to help in the design and optimisations of new OPV devices.