organic solar cells 
 ( oscs ) are among 
 the new generation photovoltaic ( pv ) technologies which could be employed 
 to convert solar energy to electricity . their main advantage over 
 more traditional pv cells is the solution processability of the layers 
 that compose the device . moreover , the potential mechanical flexibility 
 and low weight of organic pv panels may enable the introduction of 
 innovative new products . the most efficient oscs 
 require a bulk heterojunction ( bhj ) , an active layer 
 composed of intimately intermixed electron acceptors and electron 
 donors . not only their electronic properties but also the morphology of the active layer will 
 determine the final device efficiency . these two different but entangled 
 aspects of a chosen donor  acceptor system along with the many 
 degrees of freedom offered by organic materials translate into a vast 
 multiparameter optimization space which , despite extensive research 
 efforts , has not yet led to efficiencies which can compete with inorganic 
 pv technologies . design principles for bringing bhj solar cells to 
 a substantially higher performance level are still sought after . in oscs acceptor interface before 
 the resultant free carriers can be transported to the electrodes . 
 a large donor acceptor interfacial area must therefore be present 
 for efficient exciton dissociation ; second , 
 both phases need to be continuous and well - connected to the respective 
 electrode for adequate charge collection . furthermore , a high degree 
 of crystallinity is usually beneficial , as it increases charge mobility , 
 while the active layer needs to be at least 80100 
 nm thick so as to absorb an adequate amount of light . from these conditions such networks , usually achieved by spin coating 
 or doctor blading , are generally believed to be a kinetically trapped 
 state , as they are formed during the rapid solidification of a thin 
 layer of solution cast on a substrate . many processing conditions have been found critical for the resulting morphology , but determining a priori those which will lead to an optimal morphology 
 when the two constituting materials are modified is still not possible . atomic force microscopy , and transmission and scanning electron microscopy have been used to image 
 bhjs but provide only limited insights about the nature of the interpenetrating 
 network . however , the inherently 
 poor contrast due to the weak electron scattering of organic materials 
 leads to a limited and complex interpretation of the 3d - reconstructed 
 morphology . techniques used to quantify the crystallinity of bhjs 
 are grazing - incidence x - ray diffraction ( gixd ) and grazing - incidence small-/wide - angle x - ray scattering 
 ( gisaxs / giwaxs ) . the valuable information 
 provided by these techniques represents , however , averages over the 
 irradiated area of the sample , thus preventing real space localization 
 of nanocrystallites . computational modeling provides atomistic 
 detail and thereby has 
 the potential to bridge the gap between experimental and molecular 
 length scales . in particular , coarse - grain ( cg ) molecular dynamics 
 ( md ) simulations are often employed to model bhjs , being able to capture the morphology at relevant 
 length scales while still retaining a molecular description ( as opposed 
 to monte carlo simulated annealing or mean field approaches where that is lost ) . however , cg 
 models developed to date are not readily transferable . this hinders their potential use within a high - throughput screening 
 scheme aimed at optimizing the morphology of a novel donor furthermore , cg models currently applied in this field often simplify the description of donor 
 and acceptor molecules too severely : while this allows reaching length 
 scales which are comparable to real device thicknesses , the loss of detail hampers the discrimination 
 of different donor ( acceptor ) structures and does not allow for direct retrieval of atom - resolved structural information . such structures are needed as input for higher level of theory , e.g. , 
 quantum mechanical ( qm ) , calculations that address the electronic 
 states involved in the charge separation . this does not reflect 
 the morphological evolution taking place during the solution processing 
 of bhjs , hence making the link between experimental fabrication conditions 
 and simulations more elusive . importantly , lee and pao have recently 
 developed a method to generate morphologies 
 which remedy this : during a standard md simulation , the solvent is 
 slowly removed , thus simulating the experimental evaporation process . the evolution of the bhj morphology during solution processing can 
 therefore be studied , and a direct link to experimental conditions 
 can be made . in the present work , we employ cg models which 
 retain a sizable 
 degree of chemical specificity to simulate the large scale organization 
 of polymer  fullerene systems during solution processing , which 
 results in a predictive tool for bhj morphologies . the method follows 
 the one developed by lee and pao , but 
 differs from it in the nature of the cg models employed , namely finer 
 cg models based on the martini cg force field . the chemical specificity of this force field implies the 
 possibility to use it for rational design of donor ( acceptor ) compounds 
 which differ only slightly in chemical nature . moreover , martini gives 
 ready access to fully atomistic detail through backmapping , opening the way for improved qm calculations . as a first case study , we chose to simulate the poly(3-hexyl - thiophene ) 
 ( p3ht)phenyl - c61-butyric acid methyl ester ( pcbm ) 
 mixture . the validity of the method is assessed by comparing the outcome 
 of simulations to a variety of experimental data available for this 
 well - studied system . the effect of drying rate , 
 polymer molecular weight , and thermal annealing are investigated extensively , 
 resulting in trends in accordance with experimental observations . 
 finally , the atomistic detail of the whole blend is directly retrieved 
 by means of backmapping . together , the results presented in this contribution 
 show the potential of the proposed methodology to get insight into 
 larger scale morphological organization while , at the same time , obtaining 
 atom - resolved structural information on , for example , donor  acceptor 
 interfaces where such details are vital to solve the puzzle of the 
 charge separation mechanism . solvent evaporation cg md simulations 
 were performed using newly developed martini cg models of p3ht , pcbm , and chlorobenzene ( cb ) . the atomistic structures 
 along with schematic representations of the cg models are presented 
 in figure 1 . on average , 4 atoms are mapped to one cg particle , while bond , angle , and dihedral 
 potentials are used to have the cg models reproduce the structure 
 and flexibility of the corresponding atomistic models . interactions 
 have been calibrated by using experimental free energy of transfer 
 data as reference . a detailed description of the cg models , atomistic 
 models which aided the parametrization and validation of the cg models , 
 as well as the parametrization and validation procedures themselves 
 are presented in the methods section and in 
 the supporting information . all - atom structures ( left ) 
 and corresponding coarse - grain models 
 ( right ) for p3ht ( in red ; a trimer is shown ) , pcbm ( in blue ) , and 
 cb ( in green ) . the connectivity between the coarse - grain interaction 
 sites is highlighted , while the actual size of the beads is shown 
 with semitransparent spheres . in figure 2 , a prototypical 
 solvent evaporation simulation is illustrated ( see also video in the supporting information ) , and the resulting bhj 
 morphology compared to experimental data from ref ( 38 ) . in particular , figure 2b shows a comparison 
 between the morphology obtained via an evaporation simulation ( where 
 the molecular weight ( mw ) of p3ht is 8 kda ) and an energy - filtered 
 scanning electron microscopy ( efsem ) image taken by masters et al . of an as - cast p3ht  pcbm blend ( p3ht mw 
 of 54 kda ) spin - coated from cb . the experimental efsem image contains gray 
 zones which masters et al . ascribed to mixed polymer fullerene and 
 p3ht are close to each other and give rise to an average signal , which 
 thus can not be resolved . by contrast , simulations allow resolution 
 of structures up to the level of single cg particles ( the cg particle - resolved 
 morphology is shown in figure s1a ) . to 
 account for this discrepancy , and for the escape depth of secondary 
 electrons from p3ht ( found to be between around 2030  ) , a 2.5 nm thick layer taken from the 
 cg morphology has been spatially discretized in polymer , fullerene , 
 and mixed phases , following the procedure outlined in the methods section . a comparison of the resulting image 
 ( left - hand side of figure 2b ) to the experimental one ( right - hand side of figure 2b ) suggests that obtained donor 
 and acceptor domain sizes are in qualitative agreement with the experimental 
 ones . however , it is important to stress that the close agreement 
 between the experimentally and computationally obtained morphologies 
 is somewhat unexpected , given the difference in processing conditions . a first difference is in the way the solvent evaporation takes place . 
 in spin coating , the solvent evaporates from the surface of the blend , 
 whereas in our simulations solvent is taken out of the system randomly . the latter approach leads to a more uniform distribution of the molecular 
 components along the sample normal . a revised protocol in which surface 
 evaporation is taken into account is currently being developed to 
 quantify this effect in future studies . experimentally , 
 solvent evaporation takes place on the time scale of seconds , whereas 
 the simulations are necessary limited to the sub - millisecond range . the length scale of the experimental samples is also larger , with 
 typical sample areas of 1 cm , compared to 900 
 nm in the simulated samples . note that the polymer chains 
 in the simulation are smaller than those used in experiment ( 8 kda 
 vs 3050 kda ) , compensating , at least to some extent , for the 
 shorter time scale . apparently , the primary driving forces governing 
 morphology formation in real blends act on time and length scales 
 accessible by our simulations . further details about these driving 
 forces during the evaporation process are discussed below . simulation 
 of the evolution of a p3ht : pcbm bulk heterojunction 
 during solution - processing ( a ) and comparison to experimental data 
 ( b ) . the evolution of p3ht - p3ht , p3ht - pcbm , and pcbm  pcbm contacts 
 and of the solvent amount during drying are plotted ( a ) , and snapshots 
 at different times during the simulation are shown ( p3ht chains are 
 rendered in red while pcbm molecules are in blue ; solvent molecules 
 are not shown for clarity ) . the three stages of bhj morphology evolution 
 ( see text ) are indicated with different shades of gray . ( b ) visual 
 comparison between a spatially discretized morphology obtained via 
 the solvent evaporation protocol ( left - hand side ; the thin white square 
 shows the repeating unit cell ) and a spin - coated blend imaged by efsem 
 ( right - hand side ; close up of figure 7c of ref ( 38 ) , reproduced under the 
 creative commons attribution 4.0 international license ) . the evolution of the morphology during the evaporation 
 can be followed 
 by computing the number of p3ht p3ht , p3ht  pcbm , and 
 pcbm  pcbm contacts during the drying process . these have been 
 calculated as described in detail in the methods section and are plotted in figure 2a . three evaporation stages can be distinguished . within 
 the first microsecond , a first phase where molecules are still in 
 solution and thus relatively free to diffuse can be identified : during 
 this phase , a substantial increase in the number of p3ht - p3ht contacts 
 due to the aggregation of p3ht chains to form ordered nanostructures 
 ( which will be discussed more extensively in the following sections ) 
 is observed . p3ht contacts initially decrease 
 as pcbm molecules are deprived of p3ht neighboring molecules ( as these 
 aggregate ) . fullerene molecules are found to be considerably more 
 soluble in cb and do not show aggregation propensity . this is consistent 
 with the experimental solubility of p3ht in cb ( 16 mg / ml for 
 p3ht with mw of 65 kda ) being lower than the one of pcbm ( 60 
 mg / ml ) . in the second stage , which starts 
 already after 1 s when the amount of solvent is less than 60% 
 of the initial amount , the diffusion processes become much slower 
 ( time scale of hundreds of nanoseconds to microseconds ) due to the 
 increased size of p3ht nanostructures and the ever - decreasing quantity 
 of solvent , which starts to reduce dramatically the mobility of the 
 molecules . as a consequence , in this phase the number of contacts 
 between all compounds steadily increases . these aggregation processes 
 become gradually more stagnant and essentially stop when the solvent 
 weight content is about 1% . at this point , in the third and final 
 stage , the blend is in an almost frozen state where the free energy 
 barriers to even minimal rearrangement are very high . these observations 
 agree with the drying dynamics observed experimentally by gixd and 
 giwaxs / gisaxs during spin - coating or 
 printing . four drying stages are commonly individuated : 
 a first stage in which no scattering is detected and where thus the 
 molecules are considered being dispersed in solution ; a second stage 
 in which the p3ht crystals start to nucleate and grow ; a third stage , 
 also called the solvent swollen glassy state , where the mobility of the molecules is already compromised by the 
 low amount of solvent and where , therefore , no new crystals can form 
 but the existing ones continue drying ; a last stage where the film 
 is practically dry and close to the final structure . in the present 
 simulations , aggregation of p3ht is observed from the beginning , which 
 seems not to be consistent with experimental findings . however , the 
 aggregation of p3ht molecules may be happening earlier also in experiments , 
 but the p3ht crystals have to reach a certain size in order to be 
 detected by x - ray scattering techniques . moreover , the initial p3ht 
 concentration in the simulations ( 39 mg ml ) is higher than the one commonly used in experiments ( 15 
 mg ml ) , positioning p3ht already beyond the solubility 
 limit . in summary , the driving force for the ( modest ) phase 
 separation 
 observed during the solvent evaporation simulation is the self - organization 
 propensity of p3ht . this process can be hampered by a faster evaporation 
 and/or entanglement of polymer chains , as explored in the next sections . it is crucial also to realize that , due to the relatively good affinity 
 of pcbm and p3ht , large - scale phase separation does not occur . slow drying , 
 also called solvent annealing or solvent - assisted 
 annealing , has 
 been shown to allow for structure reorganization , leading the p3ht : pcbm 
 blend to approach a degree of phase separation closer to thermodynamic 
 equilibrium . lower mw ( up to 10 kda ) p3ht has been found to have a 
 stronger tendency to phase separate during the solvent evaporation 
 process than higher mw fractions , likely 
 due to an increased capability to crystallize . to investigate the 
 effect of drying rate and polymer length on the as - cast morphologies , 
 we performed two additional sets of simulations : one in which the 
 rate of solvent evaporation was varied , and another where this was 
 kept fixed while varying the polymer mw . to quantify the degree 
 of phase separation , the number of p3ht  pcbm contacts in the 
 dried morphologies has been computed . fullerene 
 contacts correlates with the extent of donor  acceptor interface 
 present in the blend . fullerene 
 contacts is obtained for a planar heterojunction morphology where 
 only one interface is present ( see figure s10 ) ; on the opposite end , the maximum amount of contacts is found in 
 a completely intermixed morphology . this has been built as a random 
 mixture , as explained in more detail in the supporting information . the numbers of contacts in the planar heterojunction 
 and intermixed extremes , respectively , have been used as references 
 to normalize the computed fraction of p3ht - pcbm contacts . using this 
 scale , figure 3a and 3b fullerene 
 contacts varies with the drying rate and with the polymer molecular 
 weight , respectively . in the latter case , normalized number of p3ht  pcbm contacts expressed in percentage 
 ( where 0 is taken as the number of contacts in a planar heterojunction 
 and 100 is the one computed for a finely intermixed morphology ) are 
 shown as blue bars . top views of representative morphologies are shown 
 for the various drying rates ( increasing drying rate from bottom left 
 to top , clockwise ) and for various polymer molecular weights ( increasing 
 polymer weight from top to bottom right , clockwise ) . faster 
 evaporation 
 does not allow for polymer crystals to grow , which therefore end up 
 having smaller dimensions ( and , as a consequence , pcbm domains remain 
 smaller as well ) . this result is consistent with the hypothesis that 
 p3ht  pcbm blends are an example of blends which do not ( liquid  liquid 
 or solid  liquid ) phase separate considerably during spin coating . even though fast drying of a p3ht : pcbm mixture 
 leads to an intimately mixed blend , simulations also predict that 
 longer polymeric chains lead to less phase - separated morphologies 
 ( figure 3b ) . slower 
 diffusion and more chain entanglement decelerate the aggregation of 
 longer polymer chains . however , the trend is not so pronounced , which 
 may be due to the limited range of mws considered in the present contribution 
 ( 1 to 8 kda ) . in fact , liu et al . reported 
 that blends with a p3ht mw up to 10 kda practically behave indistinguishably 
 during spin coating . pcbm solar cells require an annealing 
 step in order to 
 reach their full potential . the improved 
 power conversion efficiency is well documented and has been shown 
 to involve polymer crystallization  which increases optical 
 absorption and improves charge transport  and optimized phase 
 segregation , leading to more efficient charge separation . the width 
 of polymer and fullerene domains is known to increase upon annealing , 
 up to dimensions in the 1020 nm range . the mobility of pcbm molecules 
 or small aggregates has been demonstrated by several studies . it is often concluded to result from crystallization of p3ht disordered 
 domains which expel fullerene loading . in particular , diffusion of fullerenes appears to occur only through 
 the disordered regions of p3ht . conjectures 
 about polymer crystal growth during annealing inhibiting or competing with the formation of large fullerene domains are also found in the literature . 
 in order to investigate the effect of thermal annealing on the obtained 
 morphologies , a thermal annealing protocol has been applied by carrying 
 out md simulations with an increased temperature , as described in 
 detail in the methods section . figure 4a shows three snapshots 
 taken at different times during an annealing simulation for a p3ht only p3ht backbones are shown to 
 highlight the increased ordering in the p3ht phase upon annealing . the voids between p3ht sheets are , in general , occupied by p3ht side 
 chains , while the bigger white areas denote the presence of pcbm domains 
 ( compare to figure s2 ) . when p3ht is given more time for reorganization , crystallinity increases 
 and domain segregation is enhanced . the 
 simulations further provide insight into the driving force for the 
 morphology evolution upon annealing . following what is also observed 
 during the solvent evaporation simulations , p3ht drives the domain 
 segregation as p3ht crystals grow and therefore expel pcbm molecules . no breaking of p3ht stacks is observed , corroborating the hypothesis 
 that diffusion of fullerene takes place through amorphous domains . snaphots at different times during the annealing 
 process ( a ) where only the p3ht backbones are shown . white areas denote 
 the location of p3ht side chains or pcbm domains ( compare to figure s2 ) . a zoom in on the blend reveals a 
 stack of 3 p3ht chains ( b ) reporting the observed cg stacking distance 
 of 4.5  . ( c ) the corresponding computed scattering profiles 
 are shown ; gaussian fits to the ( 100 ) and ( 010 ) peaks are shown in 
 gray . scattering curves have been computed 
 as described in the methods section and are 
 shown in figure 4c . note that only thiophene 
 rings are considered in the simulated scattering curves . upon annealing , the intensities of the ( 100 ) ( or lamellar ) and ( 010 ) ( or stacking ) 
 peaks , which are fairly weak in the as - cast scattering curve , increase 
 considerably , in agreement with experimental observations . the obtained peaks can be 
 fitted with a gaussian function to obtain peak positions . the sharper 
 ( 100 ) and ( 010 ) peaks observed for the annealed blend , qualitatively 
 indicating a higher degree of crystallinity , are located at positions q  0.426  and q  1.385  , respectively . these values correspond to a lamellar distance of 14.7  and 
 a stacking distance of 4.5  ( highlighted in figure 4b ) . the discrepancy with the 
 experimental stacking distance value of 3.8  is due to the too large radius of cg beads , which does 
 not allow for a ring thickness which matches the all - atom one ; the 
 lamellar distance is instead underestimated ( 16.5  , experimentally ) . such discrepancies , inherent to the current 
 martini cg force field ( see also the methods section ) , do not , however , impact the nanomorphology evolution . the effect of p3ht mw on the annealing process has also been investigated 
 by simulated annealing of blends containing p3ht with a mw ranging 
 from 2 to 8 kda . lower mw p3ht is found to both crystallize and promote 
 a phase - separated morphology more readily , as can be seen in figure s3 , where annealing results are shown 
 for blends with varying p3ht mw ( 28 kda ) . in particular , p3ht 
 chains with mws in the 24 kda range respond similarly to annealing , 
 as quantified by the decrease in p3ht  pcbm contacts shown in figure s3 . when the mw is raised to 8 kda , a 
 lower decrease of p3ht  pcbm contacts is observed ( figure s3 ) : the blend responds less promptly 
 to thermal annealing due to a slowed down dynamics caused by the increased 
 molecular weight . moreover , in this latter case amorphous regions 
 can still be observed ( figure s3c ) , while 
 in the former case the p3ht phase is almost completely crystalline . we can therefore conclude that lower ( up to 4 kda ) mw p3ht shows a 
 higher stacking efficiency . low - mw p3ht films ( < 4 kda ) are known 
 to be more sensitive to processing conditions . those protocols which 
 give more time for molecules to approach thermodynamic equilibrium 
 ( e.g. , slow growth , thermal annealing ) lead to more crystalline phases , as the molecules diffuse relatively fast due 
 to their small size . by contrast , polymer chain entanglement and lower mobility slow the kinetics of ordering 
 and , consequently , the phase separation process in high ( > 30 kda ) 
 mw p3ht films . the present data are also 
 in agreement with gixd findings by liu et al . which show how low - mw ( up to 5 kda ) p3ht , both in pure p3ht films 
 and when blended with pcbm , shows considerable reorientation freedom 
 upon annealing , while this effect is considerably reduced already 
 at mw of 10 kda . the observed p3ht structures correspond to 
 the well - characterized 
 2d sheets in which the polymer is known to self - organize . the strong tendency of p3ht to self - organize is well - known , and it is evident from the outcome of the simulations that the model 
 also strongly favors aggregation . this self - organizing propensity 
 seems to be induced by the geometrical features of p3ht chains . it 
 should be stressed that no explicitly ad - hoc anisotropic 
 interaction has been introduced between thiophene rings , in contrast 
 to the strategy proposed by carrillo et al . in order to improve the ( also referred to as the lpc model ) . in that model , representing thiophene units 
 with single beads sacrifices one of the characteristic features of 
 this polymer , namely the planar thiophene backbone . in the present 
 model , the employed finer mapping preserves the thiophene planarity , 
 and thus anisotropy , which evidently contains the driving force for 
 p3ht self - organization . the cg morphologies obtained have finally been 
 backmapped to atomistic  also 
 termed all - atom ( aa)resolution by employing a published protocol described in the methods section . a snapshot of a cg donor  acceptor interface for 
 an annealed blend along with its atom - resolved counterpart is shown 
 in figure 5 . in this 
 procedure , the fully relaxed atomistic configuration is obtained by 
 placing atoms in positions consistent with the cg beads that represent 
 them and then relaxing the structure to the local atomistic energy 
 minimum . such a direct backmapping methodology differs 
 from the ones employed by , for example , carrillo et al . , where only a qualitative estimation of the donor  acceptor 
 interface configuration can be made from the large - scale cg simulations . 
 in that case , the obtained cg morphologies can not be unambiguously 
 backmapped to aa resolution due to too severe loss of detail inherent 
 to the lpc model . the qualitative information can only be used to 
  manually  assemble an atomistic interface . acceptor 
 interfaces , such as the one of liu et al . , where an idealized donor  acceptor 
 interface is assembled by putting together p3ht and pcbm crystals . while there are always multiple aa structures 
 corresponding to the same cg configuration , the present finer cg models 
 restrict their number and allow for the use of a semiautomatized procedure 
 for the direct retrieval of structural detail . this can pave the way 
 to qm calculations which can fully benefit from large - scale derived 
 structural information . in particular , the investigation of the local 
 dielectric response of organic photovoltaic blends upon changes in 
 the electronic density is of notable interest , as the dielectric screening 
 of charges is crucial for the exciton splitting and recombination 
 processes . promisingly , previous studies 
 on model systems have shown how the reorientation of dipole moments 
 can stabilize the charge separated state relative to the lowest charge transfer state , thereby making the 
 exciton dissociation enthalpically favorable . the coarse - grain morphology ( top ) has been backmapped to atomistic 
 resolution ( bottom ) . in both cases , only p3ht backbones ( in red ) we 
 presented a cg md based method which is able to capture the 
 kinetically trapped nature of spin - coated bulk heterojunction morphologies 
 with direct access to atom - resolved structural information . the 
 effects of slow - drying and annealing follow what is known for p3ht : pcbm 
 solar cells . the nanoscopic picture which emerges from the simulations 
 suggests that the ( moderate ) phase separation which is observed for 
 p3ht : pcbm blends is driven by the self - organization propensity of 
 p3ht . the good affinity between p3ht and pcbm does not allow for complete 
 phase separation on spin coating time scales . the resolution 
 mapping scheme of the martini cg force field represents 
 a favorable compromise between the speedup of cg models and the preservation 
 of chemical detail ; moreover , the existence of a robust backmapping 
 protocol ensures the retrieval of the atomistic structures underlying 
 the cg ones , opening the way for a molecular view on an in 
 silico solution - processed bulk heterojunction . as a 
 whole , this work proposes a predictive methodology for obtaining 
 morphologies of solution - processed soft matter systems with atomistic 
 resolution . the morphologies can subsequently serve as a starting 
 point for qm calculations on excitonic properties , or be used to compute 
 the mechanical response of bhj materials . cg md 
 simulations have been run using models 
 which are based on the martini cg force field . developed initially for biomolecules , this cg force field has been later extended to nonbiological systems , 
 including various polymers and carbon - based nanostructures . as a general rule , four heavy 
 atoms finer mapping 
 is , however , allowed for ring systems , where an interaction site may 
 represent two or three atoms . model , a 16 beads representation of c60 fullerene , while the 
 oligothiophene amphiphile martini model has been used as starting point for the p3ht cg model . experimental 
 transfer free energies were the main target of the martini force field 
 parametrization . this ( top - down ) approach 
 is used to obtain nonbonded interactions between cg particles . these 
 are determined by the type of beads which are used to describe the 
 molecules , the choice of which is made on the basis of experimental 
 free energy of transfer data . bonded parameters between the beads are instead obtained by matching bond and angle distributions of atomistic 
 simulations ( bottom - up ) . both cg and aa models have first been 
 validated by comparing computed free energies of transfer between 
 different solvents to experimental data for various molecular fragments . the results , collected in table s6 , show 
 good agreement between experiments and simulations . interactions between 
 the molecules studied here were further validated by comparing dimerization 
 potentials of mean force ( pmfs ) at the cg level with the corresponding 
 aa ones . pmfs were calculated for the dimerization of two 3ht molecules , 
 two pcbm molecules , and the 3ht - pcbm pair , all in cb solution . a thorough description of the free energy of 
 transfer and pmf calculations can be found in the supporting information . as already hinted previously , 
 the martini cg force field parametrization , 
 based on a four heavy atoms to one cg site mapping scheme , brings 
 with it a bead size which is too large to accurately reproduce the 
 thickness of ring structures . this causes the discrepancy in the cg 
 stacking distance observed in the simulated morphologies ( and also 
 seen in the potential of mean force of dimerization of two p3ht monomers 
 in the supporting information , figure s7 ) . the bead radius is part of the martini cg force field and can not be 
 tweaked arbitrarily without disrupting the carefully parametrized 
 martini partitioning equilibria . work in our research group is being 
 carried out on improving properties of martini cg ring models . nevertheless , 
 the increased thickness of thiophene rings at the cg level does not 
 affect the nature of p3ht self - assembly , which is evident from the 
 simulations and which gives rise to structures in agreement with experiments . the approach follows 
 the one recently developed by lee and pao . starting from a simulation box ( 30  30  88 nm ) containing a ternary mixture p3ht  pcbm  cb ( 39 
 mg ml in p3ht and 31 mg ml in pcbm for a 1:0.8 weight ratio ) , 1.25% of the solvent ( i.e. , cb ) 
 is removed every time interval t until a dried blend 
 is obtained ( 30  30  5 nm ) . the degrees 
 of polymerization n of ( regioregular ) p3ht employed 
 are 6 , 12 , 24 , 36 , and 48 , approximately corresponding to molecular 
 weights of 1000 , 2000 , 4000 , 6000 , and 8000 da . various sizes for 
 the system have been tested , and the simulations have been found not 
 to suffer from finite size effects ( see the supporting information ) . the largest size which could be exhaustively 
 simulated was chosen . after every solvent removal , 4 ns of npt equilibration are run followed 
 by a run 180 , 120 , 60 , 30 , 15 , 7.5 , or 3 ns long . this leads to total 
 drying times of 100 , 70 , 36 , 19 11 , 6.5 , and 4 s , respectively . the equations of motion were integrated with a time step of 30 or 
 20 fs ( the latter employed in the equilibration and often in the second 
 half of the evaporation to avoid numerical instabilities ) . the box 
 dimensions were fixed in the lateral directions by setting the compressibility 
 to 0 bar . lennard - jones potentials and forces are cut off at 1.1 nm ;  potential 
 modifiers  are used to shift the potentials to zero at the 
 cut off . the velocity - rescaling thermostat and the parrinello  rahman barostat are employed to maintain pressure and temperature , 
 respectively , with coupling parameters of 1.0 and 12.0 ( or 15.0 ) ps . all the run parameter files , 
 as well as the solvent evaporation protocol , implemented as a bash 
 script ( which assumes the use of gromacs ) , are available for download 
 as part of the supporting information and 
 on the martini portal http://cgmartini.nl . the temperature 
 has been raised gradually , in the following way : 20 ns at 498 k , 180 
 ns at 598 k , and up to 1.8 s at 698 k. this led to a total 
 simulation time of 2 s for the full annealing cycle . the parameters 
 employed for the run are the same as the ones employed in the solvent 
 evaporation simulations ( with the exception of the temperature ) . the 
 employed simulated annealing temperature is higher than the temperatures 
 commonly employed in experiments ( 420 k ) . however , time scales are also different , as blends are 
 usually annealed for time scales currently not available for cg modeling 
 ( 510 min ) . the distances between 
 all the n thiophene centers of mass have been computed . the n by n obtained distance matrix 
 is , evidently , symmetric and thus contains n(n  1)/2 unique distances . these distances have been 
 binned in a histogram ( occurrence vs distance ; bin width = 1 nm ) . this is done in the following way : the atoms are considered 
 as point scatterers located at points in the real space given by the 
 coordinates obtained by the md simulations . if so , the scattering 
 intensity i(q ) in the 3d reciprocal 
 space is given by1where q is the reciprocal space 
 vector , rj is the position vector 
 of atom j , and zj is the atomic number of atom j. in this particular 
 case , we will make use of a one - dimensional version of eq 1 . moreover , as we consider identical 
 atoms ( i.e. , the thiophene centers of mass ) , the factor zj can be omitted . the procedure , which 
 makes use of the mdanalysis package , has been implemented 
 in a jupyter / ipython notebook available for download as part of the supporting information . briefly , mapping files , i.e. , initial positioning of atoms expressed 
 on the cg particles space , have been defined for p3ht and pcbm . the 
 program backward.py places the atoms according to 
 the definitions contained in the mapping files . a series of energy 
 minimizations is then carried out ( as implemented in the bash script initram.sh ) until a relaxed atomistic morphology is obtained . the programs are available for download as part 
 of the supporting information and on the 
 martini portal http://cgmartini.nl . the number of contacts 
 between molecules a and b in the simulations were calculated by counting 
 one contact for each cg bead belonging to a molecule b within a sphere 
 of radius 0.6 nm from each cg bead belonging to a molecule a. the 
 cutoff distance of 0.6 nm comprises the nearest neighbor cg sites 
 around a cg particle , with the radius of martini s - particles being 
 0.24 nm ( where n = 0.43 nm ) . when counting contacts , the 
 counting of intramolecular 
 contacts is avoided . more details on the computation of the number 
 of contacts and generation of the planar heterojunction and intermixed 
 morphology used to normalize the number of contacts can be found in 
 the supporting information . from each obtained dried bhj ( of dimensions 30  30  5 nm ) , two slabs ( of dimensions 30  30  2.5 
 nm ) have been extracted . the thickness of these slabs 
 ( 2.5 nm ) , and consequently of the voxels into which the slabs are 
 divided , has been chosen based on the p3ht escape depth of secondary 
 electrons in the efsem measurements ( found to be between around 2030 
  ) to which the spatially discretized 
 morphologies are compared . the morphology slabs have been therefore 
 converted to a grid of polymer ( red ) , fullerene ( blue ) , and mixed 
 ( gray ) domains employing the following algorithm : the simulation box 
 is divided into voxels of dimensions 1  1  2.5 nm ; the number of p3ht ( n ) and pcbm ( n ) particles is then computed for each i , j voxel , and a color assigned to the 
 voxel based on the relative number of polymer and fullerene cg particles 
 contained in the voxel . in particular , the following ratio is evaluated : a xi , j fraction 
 higher ( lower ) than 0.6 ( 0.6 ) 
 has been assigned to the polymer ( fullerene ) phase , while fractions 
 between 0.6 and 0.6 ( included ) have been assigned to the mixed 
 phase .
control 
 over the morphology of the active layer of bulk heterojunction 
 ( bhj ) organic solar cells is paramount to achieve high - efficiency 
 devices . 
 however , no method currently available can predict morphologies 
 for a novel donor 
 acceptor blend . 
 an approach which allows 
 reaching relevant length scales , retaining chemical specificity , and 
 mimicking experimental fabrication conditions , and which is suited 
 for high - throughput schemes has been proven challenging to find . here 
 , 
 we propose a method to generate atom - resolved morphologies of bhjs 
 which conforms to these requirements . 
 coarse - grain ( cg ) molecular 
 dynamics simulations are employed to simulate the large - scale morphological 
 organization during solution - processing . 
 the use of cg models which 
 retain chemical specificity translates into a direct path to the rational 
 design of donor and acceptor compounds which differ only slightly 
 in chemical nature . finally , the direct retrieval of fully atomistic 
 detail is possible through backmapping , opening the way for improved 
 quantum mechanical calculations addressing the charge separation mechanism . 
 the method is illustrated for the poly(3-hexyl - thiophene ) ( p3ht)phenyl - c61-butyric 
 acid methyl ester ( pcbm ) mixture , and found to predict morphologies 
 in agreement with experimental data . 
 the effect of drying rate , p3ht 
 molecular weight , and thermal annealing are investigated extensively , 
 resulting in trends mimicking experimental findings . 
 the proposed 
 methodology can help reduce the parameter space which has to be explored 
 before obtaining optimal morphologies not only for bhj solar cells 
 but also for any other solution - processed soft matter device .