Daniele
Padula
*a,
Alessandro
Landi
*b and
Giacomo
Prampolini
*c
aDipartimento di Biotecnologie, Chimica e Farmacia, Università di Siena, Via A. Moro 2, 53100 Siena, Italy. E-mail: daniele.padula@unisi.it
bDipartimento di Chimica e Biologia Adolfo Zambelli, Università di Salerno, Via Giovanni Paolo II, I-84084 Fisciano (SA), Italy. E-mail: alelandi1@unisa.it
cIstituto di Chimica dei Composti OrganoMetallici (ICCOM-CNR), Area della Ricerca, Via G. Moruzzi 1, 56124 Pisa, Italy. E-mail: giacomo.prampolini@pi.iccom.cnr.it
First published on 20th July 2023
Recent advances in non-fullerene acceptors (NFAs) have significantly increased the efficiency of organic photovoltaics, reaching approximately 20% in single junction solar cells. These advancements are attributed to the introduction of the promising L8-R series, derived from the high-performance NFA known as Y6. To gain a deeper insight about such efficiency improvement, here we computationally characterise Y6 and its derivatives, focusing on the supramolecular structure of the resulting aggregates and their electron transfer properties. The applied computational protocol indicates an evident relationship between different side chains and the electronic features of the NFA supramolecular architectures, and provides a possible rationale for their impact on photovoltaic device efficiency. Using density functional theory and experimental crystal structures, we show that shorter and branched alkyl side chains on L8-R derivatives improve electron transfer integrals and charge mobilities compared to Y6. This improvement is attributed to differences in crystal packing. The effect of thermal fluctuations on the performance of the NFA aggregate is instead investigated through molecular dynamics simulations with quantum-mechanically derived force fields. Results indicate that, despite a minimal impact on electron transport capabilities due to dynamic disorder, the various substitution patterns significantly influence the supramolecular arrangement of the aromatic cores. The validation of the presented computational protocol, which integrates in a multi-level procedure accurate charge transfer rates computed on reliable morphologies obtained through classical molecular dynamics with refined force fields, paves the way towards a bottom-up modelling of donor/acceptor interfaces with polymeric donors, where the accurate description of structural and electronic disorder is key to reach a computationally-driven identification of higher performing components for organic photovoltaics.
The major challenge in OPVs is achieving high power conversion efficiency (PCE, η), a critical point for commercialisation.1 After 2012, the efficiency of OPV cells reached a plateau around 10–12%, a level of performance unchallenged until very recently,2,4,5 when researchers introduced a small-molecule non-fullerene acceptor (NFA) called Y6, which enabled OPV cells to achieve efficiencies exceeding 15%.6 Since then, Y6 has become a key component of many high-efficiency OPV cells,7–13 typically in combination with polymeric donors but also with small molecules.14 Following its discovery, Y6 garnered significant attention and both theoretical and experimental studies have been devoted to rationalise and even further improve its exceptional features.15–21
From a theoretical point of view, researchers traced back Y6 excellent charge and exciton transport properties to its strong delocalisation capability, which results in a reduced chance of recombination.15–17 On the experimental side, while it is debated whether Y6 alone could lead to photocurrent and hence function as a photovoltaic material without a heterojunction,21 several derivatives were proposed, including a systematic study for a set of NFAs with branched alkyl side-chains (see Fig. 1) employed in heterojunctions with a polymeric donor (PM6).18 Indeed, chemical substitution of different side chains decorating the Y6 core led to unparalleled efficiencies in solar cells, with better performances registered for the Y6 derivative with shorter branched side chains, which showed efficiencies up to 18.3%.20
Fig. 1 Structure of small-molecule acceptors with varying side-chains and moieties constituting the A–D–A′–D–A architecture highlighted. |
Such encouraging results seem to suggest that different crystal packing induced by the modification of Y6's alkyl pendants might be the cause of such significant improvement of charge and exciton transport. As a matter of fact, differences in alkyl side chains (of both acceptor and donor components) are known to impact the performance of the resulting OPV devices,22–29 being pivotal to modulate both cristallinity30 and solubility in organic solvents.31 This, in turn, is a central point to improve the processability of the device. For instance, beside enabling solution printing,32 it affects the miscibility between components, increasing the chance of charge transfer from the donor to the acceptor, without amplifying the recombination processes that could alter charge migration.33,34
To gain deeper insight onto the mechanisms connecting chemical substitution patterns with device performances, in this work we study the set of L8-R NFAs experimentally investigated in ref. 20, focusing on the possible differences induced by the varied alkyl chains on the crystal packing and its consequent electronic behaviour. In this respect, it should be stressed that such differences are indirect effects, because alkyl side chains are not expected to affect the electronic structure of the π-conjugated core to a significant extent. The ability to predict such features is in fact critical to shift the paradigm from a trial-and-error, experimentally-driven research to a rational design based on theoretical predictions.35 Yet, the effect of side chains on OPV performance is not easily predictable: although one can speculate that in similar supra-molecular architectures π-conjugated cores substituted with longer alkyl chains possibly experience stronger effects from an increased number of low frequency vibrations,36 there are cases where side chains of intermediate length resulted in more efficient devices,29 or where longer side chains induce extended planarity in the π-conjugated core, resulting in a performance increase.37
To accomplish the challenging task of an in silico design of efficiently substituted NFAs, three main goals need to be achieved: (i) a prediction of morphology (or crystal structure) of the components38–42 (ii) a prediction of electronic properties of the aggregate in its morphology (iii) a prediction of the morphological and electronic interaction between components. In previous work, some of us have shown how point (iii) could be approached studying donor/acceptor blends based on small molecules.17 Here we instead focus on point (ii), and apply our computational protocols to the set of L8-R NFAs experimentally employed in ref. 20. By investigating their aggregation and morphology, we aim at giving a rationale on the impact of the resulting supramolecular structures on the predicted charge and exciton transport properties. In particular, on the one hand we will exploit the specific chemical structure of the considered NFAs (which all share the same core but differ for their alkyl pendants, see Fig. 1) and the homogeneous composition of the simulated NFA aggregates, to simplify the comparison among the considered compounds and clarify the relationship between the different structures and the computed electronic properties. On the other hand, it should be noticed that a positive performance of the computational protocols here proposed for simpler homogeneous systems would pave the way to future investigations on binary or ternary blends, as for instance explicitly including polymeric donors.
All four structures were employed to compute the transfer integrals between frontier orbitals as
Jij = 〈ϕi||ϕj〉, | (1) |
Together with additional geometry optimisation and frequency calculations (which we carried out on the π-conjugated core only, assuming a similar behaviour for the four NFAs), transfer integral values allow to evaluate electron transfer rates according to Fermi's Golden Rule (FGR),
(2) |
Unless otherwise stated, all electronic calculations were carried out with the Gaussian16 software,56 within the Density Functional Theory (DFT) framework, making use of the B3LYP functional and the popular 6-31G* basis set, accounting for dispersion through the Grimme D3 correction scheme.57 We chose FGR to evaluate the rates because of its improvements with respect to the widely used classical Marcus theory, since it allows taking into account tunneling effects and frequency changes between initial and final electronic states, including the whole sets of normal modes in the computation,52,53,55,58,59 leading to extremely reliable rates of radiative and non-radiative electronic transitions.54,59–61
Predicted electron transfer rates were used to estimate the drift mobility of electrons μ according to Einstein–Smoluchowski's relation
(3) |
(4) |
It should be noted that eqn (4) implies the adoption of a hopping transport mechanism, whose applicability to organic crystals where J ≈ λ is highly debated.62–66 However, both theoretical estimations67 and experimental measurements68 of mobility carried out at different temperatures show that the hopping model is adequate for the evaluation of charge mobility, at least at room temperature, in particular when quantum mechanical effects are included in the rates,58,59,69 such as in the case of the FGR used here.
Additionally, notwithstanding other popular protocols employed in the prediction of charge mobilities, such as the transient localisation theory, are generally valid for any molecular crystal, they have been usually used for perfect crystals showing herringbone structures.70 Finally, we also highlight that eqn (4) is valid only for isotropic systems,71 thus it is not entirely compatible with classical MD simulations (vide infra).
The first strategy has the clear advantage of precisely identifying the lattice or molecular vibrations directly responsible for large modulations of transfer integrals, but as a downside it assumes the harmonic approximation, which is valid only for small displacements from the equilibrium geometry. This is an unlikely event for low frequency modes, which usually show a rather flat potential energy profile, often characterised by several local minima. The second strategy instead relies on retrieving the large modulations of transfer integrals from the spectral density obtained by Fourier Transforming the autocorrelation function of the time series to a frequency domain description. As such, no direct assignment of vibrations is possible, although anharmonicity is accounted for.
While in future work we plan to compare the two strategies, in this work we resort to the classical MD approach. The main drawback of MD based protocols is the possibly poor description of the potential energy surface provided by general purpose force fields (FFs). A possible route to overcome this lack consists in resorting to highly accurate quantum mechanically derived force fields (QMD-FFs), which can be specifically tailored for the molecules under study. Indeed, the benefits of adopting such accurate description have been recently demonstrated in the domain of spectroscopic and electronic properties of organic semiconductors.81,82 Furthermore, the importance of the QMD-FF specialised description for organic aggregates was also recently assessed through simulations of the self-assembly of similar small organic molecules into orientationally ordered liquid crystal phases,83–85 where the phase transition is modulated by a delicate interplay of intramolecular flexibility and intermolecular interactions.
In fact, as detailed in Section S2 in the ESI,† the OPLS general-purpose FF does not reproduce well the quantum chemical interaction potential energy surface computed at DFT level, yielding an average error larger than 10 kJ mol−1. Even limiting to the most probable configurations, hence comparing selected interaction energy profiles (see Fig. S3 and S4, ESI†), minima are significantly shifted, and the relative stability of different molecular arrangements is not correctly described.
This poor description of intermolecular interactions could in principle significantly affect packing, and consequently transfer integrals. For these reason, we first fitted all atomic point charges, according to the RESP procedure,92 to the B3LYP-D3/6-31G* electrostatic potential of the entire molecule, and thereafter parameterised the intermolecular Lennard-Jones terms of the aromatic core following the Picky procedure as detailed in the ESI.† In summary, the QMD-FF parameterisation of the intra- and inter-molecular terms were respectively carried out with the Joyce93 and Picky94 packages. Additionally, to separately analyse the impact of point (i) and (ii) to the description of dynamic disorder through MD simulation, two other FFs were assembled for each considered molecule, according to the procedures summarised in Table 1 and detailed in Section S1.2 of the ESI.†
Fragment | Parameters | ||
---|---|---|---|
π intramolecular | OPLS | Joyce | Joyce |
π intermolecular | OPLS | OPLS | Picky |
Alkyl intramolecular | OPLS | Joyce | Joyce |
Alkyl intermolecular | OPLS | OPLS | OPLS |
Label | OPLS | Hyb-FF | QMD-FF |
Fig. 2 Unique charge transfer pathways and associated transfer integrals J identified in the crystal structure of L8-BO (L8-HD and L8-OD displays the same pathways). Side chains omitted for clarity. |
The Y6 crystal is instead different from the ones of the L8 series, although the same three charge transfer pathways can be identified. Table 2 reports a characterisation of the charge transfer pathways computed for the investigated crystals, resulting in a good agreement with electron transfer integrals and rates computed on crystal geometries and aggregates of Y6 by others.15,96 From the data reported in Table 2 we already see a trend of electron mobility that parallels the one of PCE in the corresponding heterojunctions with PM6: shorter and branched side chains lead to higher electron mobilities. The difference in the mobility between the species under study is easily rationalisable by looking at the transfer integral values of the interacting pairs. Indeed, we can see that J2 and J3, the transfer integrals for the pairs with the highest distance r (i.e. the ones which most influence the diffusion coefficient) are significantly higher for L8-BO and this has a trivial impact on the final mobility.
NFA | Path | r/Å | J/meV | k/s−1 | μ/cm2 V−1 s−1 |
---|---|---|---|---|---|
Y6 | J 1 | 9.4 | −59 | 3.75 × 1013 | 1.88 |
J 2 | 18.5 | 18 | 3.49 × 1012 | ||
J 3 | 14.0 | −15 | 2.46 × 1012 | ||
L8-BO | J 1 | 13.8 | −4 | 1.70 × 1011 | 8.11 |
J 2 | 18.8 | 45 | 2.23 × 1013 | ||
J 3 | 13.8 | −82 | 7.32 × 1013 | ||
L8-HD | J 1 | 12.2 | 5 | 3.06 × 1011 | 2.13 |
J 2 | 18.9 | 22 | 5.01 × 1012 | ||
J 3 | 13.7 | −43 | 1.98 × 1013 | ||
L8-OD | J 1 | 12.7 | −18 | 3.39 × 1012 | 2.02 |
J 2 | 18.5 | 17 | 3.06 × 1012 | ||
J 3 | 12.8 | 46 | 2.29 × 1013 |
These results seem to indicate that one of the reasons of the better performance of L8-BO with respect to the related NFAs can be the different transport properties in the microcrystalline acceptor domains that form in the donor/acceptor blend in the solar cell. Of course, this hypothesis has to be confirmed through accurate studies of the blends, a task we are currently working on and that we plan to report in forthcoming work. This explanation is in agreement with qualitative analyses based on visual inspection of the L8-R crystal structures, suggesting L8-BO superiority due to multiple charge transport channels, as well as with experimental data of electron mobilities in the blend.2,20 However, it is well known that transfer integrals are not the only parameter affecting transport. In fact, as discussed in the previous section, a crucial role is also played by dynamic disorder,72 whose main contribution is given by electron–phonon coupling. To evaluate this effect and the relative role of transfer integrals and dynamic disorder in ultimately determining the experimentally observed trend, we resorted to classical MD simulations.
In Table 3 we report the results obtained with MD simulations carried out at 150 K, i.e. the temperature at which the experimental crystal structure characterisation was performed. Further computational characterisation at 100 K and 293 K is available in the ESI.† On the one hand, as far as the structural parameters of the unit cell are concerned, crystallographic axes and angles seem not to vary much and negligible differences arise among the benchmarked FFs, all three yielding values in rather good agreement with the experimental measures, with errors ranging from less than 1% to ≈7%. On the other hand, density and volume are known to be observables more sensitive to intermolecular interactions, hence reflecting in the better agreement with experimental values registered for a complete reparameterisation (QMD-FF). Interestingly, partially refining the general purpose (OPLS) FF in its intra-molecular term (Hyb-FF) also leads toward a more realistic density. This further confirms that the overall description of the aggregate is regulated by a delicate interplay between internal flexibility and intermolecular interactions. Finally, the order of the crystal expressed in terms of the second rank order parameter P2 directed either along the largest or the smallest eigenvector of the Saupe ordering matrix (see Section S3 in the ESI†) is comparable with all FFs employed.
Property | Expt.20 | QMD-FF | Hyb-FF | OPLS |
---|---|---|---|---|
a/Å | 27.704 | 28.46 ± 0.05 | 28.04 ± 0.03 | 27.30 ± 0.08 |
b/Å | 20.855 | 21.05 ± 0.02 | 21.43 ± 0.03 | 20.50 ± 0.03 |
c/Å | 28.363 | 26.24 ± 0.07 | 26.54 ± 0.01 | 28.8 ± 0.1 |
α/° | 90.0 | 90.0 | 90.0 | 90.0 |
β/° | 105.949 | 106.95 ± 0.02 | 106.584 ± 0.009 | 104.99 ± 0.06 |
γ/° | 90.0 | 90.0 | 90.0 | 90.0 |
V/Å3 | 15756 | 15721 ± 8 | 15950 ± 13 | 16149 ± 13 |
ρ/kg m−3 | 1248 | 1254 ± 1 | 1235 ± 1 | 1217 ± 1 |
P ‖core2 | 0.942 ± 0.001 | 0.942 ± 0.001 | 0.952 ± 0.002 | |
P ‖core4 | 0.813 ± 0.003 | 0.815 ± 0.003 | 0.845 ± 0.005 | |
P ⊥core2 | −0.500 | −0.500 | −0.495 | |
P ⊥core4 | 0.374 | 0.374 | 0.364 ± 0.001 |
Another popular way to analyse supramolecular architectures is through radial distribution functions (RDFs). In our case we partitioned the π-conjugated core in three subsystems constituted by the central D–A′–D moiety (labelled as C), and the side A moieties (labelled as S). We used this partition to compute RDFs between the centres of mass of any pair of C or S moieties. These results are reported for MD trajectories at 150 K in the central panel of Fig. 3. While in the RDFs obtained with either QMD-FF and Hyb-FF we observe essentially the same comb-like profile, typical of highly ordered aggregates such as crystals, for OPLS we obtain a much broader and less structured profile, suggesting the formation of a less ordered aggregate that can be easily disrupted.
In particular, we highlight how each peak in the RDFs obtained with refined FFs is associated to a specific intermolecular interaction: we assign the first peak at r ≈ 4 Å to π-stacking interactions between S moieties belonging to different molecules (those seen in J1 and J3 couplings, see Fig. 2), while the second peak at r ≈ 6 Å is due to interactions between S and C moieties (those in J2 couplings, see Fig. 2), and all other peaks can be assigned easily. Conversely, in the RDF obtained with the general purpose FF, only the first peak is assignable, while at longer distances there are several overlapping contributions, diagnostic of a less structured and stable aggregate. To confirm this finding, the thermodynamic stability of the resulting NFA aggregates was first investigated through comparison of crystal parameters obtained from simulations at room temperature (the operational condition of solar cells), whose results are reported in Table S1 in the ESI.† Thereafter, heating simulation sequences (see Section S1.3.1 in the ESI†) were carried out in the 100–600 K range, and the resulting RDFs at different temperatures are displayed in the bottom panel of Fig. 3.
We notice that the general purpose description starts deviating significantly from the ordered behaviour of a crystal at an early stage of the heating sequence: the peak at r ≈ 4 Å is significantly reduced already at about room temperature, and completely lost above T ≈ 400 K, providing RDFs reminiscent of isotropic aggregates. On the contrary, the RDF peak broadening, induced by the increasing, temperature-driven disorder, is significantly slower for the QMD-FF and longer simulations are needed to observe a net transition to an isotropic phase.
Although a quantitative determination of transition temperatures would require much more extensive simulations,83,85 beyond the scope of the present work, as a general trend it appears that either partial (Hyb-FF) or total (QMD-FF) reparameterisation allow for a more accurate and reliable description of the crystal.79,100 More importantly, a general purpose representation provides a description of the aggregate in poor agreement with real-world experiments, both quantitatively (with an estimated L8-BO melting temperature TOPLSm < 530 K that badly compares with experiments yielding Tm ≈ 593 K) and qualitatively, since the crystal completely loses its order within few nanoseconds (see Fig. S11 in the ESI†), upon heating, while phase transitions for organic molecules of comparable size require significantly longer simulations.83,85
In conclusion, all the above discussed validation tests agree in indicating that at least a partial refinement of general purpose descriptions is required to achieve reliable results.79,100 Besides the intermolecular parameterisation, which directly impacts the relative arrangements of neighbouring NFA units and hence the charge transport properties of the resulting aggregate, it should be stressed that important differences were found also as a consequence of the sole intramolecular re-parameterisation. In this respect, the effect of the refinement of the QMD intra-molecular parameters is twofold. On the one hand, the accurate fitting of the flexible chain dihedrals rules the conformational dynamics within the chain, in turn affecting the supramolecular structure. On the other hand, a reliable representation of valence bonds and angles, based on higher level QM data, is beneficial to the successive electronic calculations carried out along the classical trajectories, allowing to overcome the so-called “potential energy surface mismatch” problem.101–103
NFA | J 1 ± σ/meV | J 2 ± σ/meV | J 3 ± σ/meV |
---|---|---|---|
Y6 | −42 ± 30 | 34 ± 20 | 0 ± 0 |
L8-BO | −7 ± 15 | −78 ± 25 | 50 ± 19 |
L8-HD | 25 ± 21 | −58 ± 31 | 15 ± 16 |
L8-OD | 35 ± 36 | −44 ± 24 | 7 ± 12 |
Such marked variation with respect to the static crystals can be traced back mainly to two factors (see Fig. 4):104 a change in the intermolecular π distances between the interacting molecules due to a displacement perpendicular to the aromatic planes (Fig. 4, bottom left panel) and a relative slide between the molecular planes (Fig. 4, upper panels) or rotation (Fig. 4, bottom right panel) around the plan normal axis.
Fig. 4 Lateral (red arrows) and vertical (blue arrow) relative displacements of two molecular planes. In magenta the relative rotation of the molecular planes. |
Before addressing a deeper discussion on the role played by such kind of supramolecular displacements on dynamic disorder, it is appropriate to again verify the FF reliability in their description. Fig. S3 and S4 in the ESI† show the remarkably increased agreement with respect to higher level DFT descriptions gained with the QMD re-parameterisation of the intermolecular FF terms. Interestingly, such test also reveals different shapes for the FF interaction energy curves computed varying either the stacking distance or the angle between the aromatic planes, the former being characterised by a rather steep and well defined minimum, while the latter presents a more flat profile. This in turn reflects on the distributions of dimer arrangements populated along the MD trajectories, where the intermolecular distance is essentially unaltered with respect to the value found in the crystal, hence suggesting that the main cause for the different transfer integrals values can be related to a relative slide or rotation, which are favourable to maximise π-stacking interactions without altering the distance between the aromatic planes.
Indeed, for such large conjugated molecules, the molecular orbitals (which determine the transfer integral values, see eqn (1)) show a complex nodal structure (see Fig. S18 in the ESI†) so that even a slight shift in their relative position and overlap can have a strong impact on the interaction among the aromatic cores. In this respect it should be stressed that, for a reliable sampling of such displaced dimer arrangements, it is of the foremost importance to rely on refined QMD-FFs, able to entail the QM accuracy in the description of complex π–π interaction energy surfaces.91,105
The difference in the transfer integrals with respect to the perfect crystal case is expected to translate in different charge mobility values. As a first estimation, we have computed the diffusion coefficient following eqn (4), by considering all the neighbouring couples in the MD system. We remark that this approach, although not rigorous because it is applicable only to isotropic systems, represents an estimation of the results from a stochastic charge hopping simulation (e.g. via kinetic Monte Carlo),3,71 and it is accurate enough to allow us to rationalise the observed experimental trend.
Indeed, inspection of Table 5 shows that the average mobility evaluated along the MD trajectory is slightly higher than the one evaluated for the perfect crystal structure. Nevertheless, the trend (L8-BO > L8-HD > L8-OD > Y6) is preserved. This trend seems to suggest that L8-R derivatives have better transport properties than the parent material and that shorter and branched side chains lead to enhanced charge mobility. As a consequence, devices built exploiting such derivatives are expected to yield enhanced efficiencies. After optimising the NFA, researchers further investigated the effect of blend morphology by studying ternary solar cells,106 where a new auxiliary donor designed based on the Y6 core was added.107 This led to unprecedented efficiencies of ≈19.3% in single-junction cells. However, it is not guaranteed that separate optimisation of the components can lead to the best results, as their synergy must also be considered.108,109
Y6 | L8-BO | L8-HD | L8-OD | |
---|---|---|---|---|
Crystal | 1.88 | 8.11 | 2.13 | 2.02 |
MD | 3.92 | 9.18 | 5.60 | 4.67 |
Concerning the latter calculations, the thorough set of MD runs and analyses of the resulting trajectories indicates that a reparameterisation of force fields for classical molecular dynamics based on quantum chemical reference data leads to much more realistic description of aggregates of organic semiconductors with respect to general purpose parameters.
Turning to the predicted properties, charge mobilities computed in single crystal and along the MD trajectory show that, in agreement with previously reported experiments, L8-R NFAs show enhanced electron transport capabilities with respect to the parent Y6 compound, with a clear trend showing that shorter and branched chains lead to the best performance. This is due to an increase in electron transfer integrals, which results from the altered crystal packing induced by the alkyl side chains.
This work is a necessary first step towards studying the performance of the NFA under study in binary or ternary blends with polymeric donors,20,106 where the importance of structural disorder on the electronic properties of these substrates110 calls for the adoption of accurate descriptions of both structural flexibility of all involved species, and supramolecular interaction patterns among all components of the blend. The good agreement between our analyses and the available experimental data for the investigated aggregates strongly supports the extension of the present protocols to multi-component systems, a work currently in progress in our group.
Footnote |
† Electronic supplementary information (ESI) available: Additional computational details, interaction energy curves, further structural and electronic characterisation along MD simulations. See DOI: https://doi.org/10.1039/d3ya00149k |
This journal is © The Royal Society of Chemistry 2023 |