Organic photovoltaics: determining the impact of the 1-chloronaphthalene additive in the PM6:Y6 blends from multi-scale simulations

Abstract

The development of PM6:Y6 active layer has contributed to significant advances in the field of organic photovoltaics. During its fabrication, using solvent additive such as 1-chloronaphthalene (1-CN) at a specific concentration has often been found to play a critical, positive role. There is increasing consensus that 1-CN molecules can tune the global morphology of PM6:Y6 active layer and thus device efficiency by slowing down the solvent evaporation process. However, it remains poorly understood how these 1-CN molecules impact the nano-scale molecular packing and electronic properties of PM6:Y6:1-CN blends (aspects that are challenging to characterize experimentally) as a function of 1-CN concentration in solution. Changes in “global” morphology correspond to modifications in size, connectivity, and crystallinity of PM6 and Y6 domains while changes in “local” morphology correspond to modifications in intermolecular interactions among PM6, Y6, and 1-CN moieties. The latter in turn impact the electronic properties related to hole/electron transfer rates between PM6/Y6 and 1-CN, electron/hole transfer rates between adjacent Y6/PM6 molecules/chains, energetic distributions of interfacial charge-transfer (CT) electronic states, and non-radiative recombination rates and corresponding voltage losses from the CT states to the ground state. By thoroughly investigating these aspects in PM6:Y6 blends with different 1-CN concentrations, via combining density functional theory calculations, all-atom molecular dynamics (MD) simulations, and coarse-grained MD simulations, we are able to rationalize how 1-CN improves the photovoltaic parameters and thus efficiencies of PM6:Y6-based solar cells, and to point out the requirements in terms of electronic properties of solvent additives for further improvement in device efficiency.

---

Broader context

The development of PM6:Y6 active layer has re-energized the community of organic photovoltaics. Especially, the structural modifications of Y6 and/or PM6 together with the morphological optimizations of corresponding active layers have significantly improved the device power conversion efficiencies over 20%. During the fabrication of these state-of-the-art active layers, solvent additives have often been used at a specific concentration and found to play a critical, positive role. While there is increasing consensus that additive molecules can tune the global morphology of active layers and thus device efficiency by slowing down the solvent evaporation process, it remains poorly understood how they impact the nanoscale molecular packing and electronic properties of active layers, since accessing them remains experimentally challenging as we are dealing essentially with nanoscale morphology at interfaces. Solving this issue calls for the recent development of robust computational methodologies that combine long-range corrected density functional theory calculations, all-atom molecular dynamics simulations, and coarse-grained molecular dynamics simulations. Using this multi-scale procedure, we drew a comprehensive picture that describe how the morphological and electronic properties of PM6:Y6 blends evolve as a function of 1-chloronaphthalene concentration in solution, and pointed out the requirements regarding electronic properties of solvent additives for further efficiency improvement.

Introduction

Organic solar cells (OSCs), as a class of photovoltaic devices that directly convert sunlight into electricity, have attracted considerable attention owing to a number of advantages including, for instance, light weight, short payback time, flexibility, solution processibility, and semitransparency.^1–3 Since the development of the Y6 series of non-fullerene small-molecule acceptors,⁴ the OSC power conversion efficiencies (PCEs) have significantly improved.^5–31 Advances have commonly involved a number of strategies, such as structural modifications of the Y6 acceptor^{9,15–17,22,23,25,27} and the polymer donors,^9,10,21 tuning of the active-layer deposition techniques,^13,14,23,30 incorporation of ternary components,^9,17,23,27 and/or introduction of solvent additives.^{6–8,10–14,16,18–20,22–24,26–31}

Adding solvent additives represents a low-cost and, as it turns out, efficient way of PCE improvement. A widely used additive is 1-chloronaphthalene (1-CN).^{6–8,10,11,13,16,19,20,23,24,31} In particular, the use of 1-CN has recently contributed to achieving 20.8% PCEs for layer-by-layer fabricated OSCs.¹³ It has been extensively reported that introducing 1-CN into the processing solution at a specific concentration (e.g., 0.5%, v/v) usually increases the short-circuit current density (J_SC) and fill factor (FF), somewhat decreases the open-circuit voltage (V_OC), and overall increases the PCEs of OSCs based on Y6 or its derivatives.^{4,11,19,20,24} Similar results are often found as well when using other solvent additives in recent state-of-the-art OSCs.^{12,18,22,26,29}

Experimental and theoretical efforts have pointed out that: (i) its high boiling point and low vapor pressure provide 1-CN with the ability to extend the solvent-evaporation process; this then gives the Y6 molecules enough time to rearrange and adjust their packing and interactions with the chains of the polymer donor paired with it (for instance, PM6, which is widely used in conjunction with Y6);^8,10,11 (ii) adding 1-CN is able to tune the sizes of the Y6 domains;^6,7 (iii) the π–π interactions between 1-CN and Y6 molecules can enhance the intermolecular interaction strength among adjacent Y6 molecules, which reduces the Y6–Y6 packing distance and improves the Y6–Y6 packing order;^6–8,10,31 and (iv) because of its non-volatile characteristics, it is difficult to remove completely the 1-CN molecules after the evaporation of the solvent molecules, which leads to a residual presence of 1-CN in the resulting active layer.^31–33

However, a number of aspects remain not well documented, in particular (i) how do the morphology of the polymer:Y6 blends at the nanoscale (“local” level) and the related electronic properties evolve as a function of 1-CN concentration in solution; and (ii) whether the residual 1-CN molecules participate in the charge transport processes taking place in the active layer. Drawing a comprehensive picture of these points would significantly improve our understanding of the important role played by 1-CN and guide further design of more efficient solvent additives. Here, we do so by taking the PM6:Y6:1-CN blends as representative systems (see Fig. 1) and by thoroughly investigating them via multi-scale simulations that combine long-range corrected density functional theory (DFT) calculations, all-atom molecular dynamics (AA-MD) simulations, and coarse-grained molecular dynamics (CG-MD) simulations.

Fig. 1 Chemical structures of poly[(2,6-(4,8-bis(5-(2-ethylhexyl-3-fluoro)thiophen-2-yl)-benzo[1,2-b:4,5-b′]dithiophene))-alt-(5,5-(1′,3′-di-2-thienyl-5′,7′-bis(2-ethylhexyl)benzo[1′,2′-c:4′,5′-c′]dithiophene-4,8-dione))] (PM6), 2,2′-[[12,13-bis(2-ethylhexyl)-12,13-dihydro-3,9-diundecylbisthieno[2′,3′:4′,5′]thieno[2′,3′:4,5]pyrrolo[3,2-e:2′,3′-g][2,1,3]benzothiadiazole-2,10-diyl]bis[methylidyne(5,6-difluoro-3-oxo-1H-indene-2,1(3H)-diylidene)]]bis[propanedinitrile] (Y6), and 1-chloronaphthalene (1-CN).

Long-range corrected DFT calculations are used to parametrize the AA-MD force field and examine the electronic properties related to the local morphology; AA-MD simulations, to parametrize the CG-MD force field and describe the local molecular packing patterns; and CG-MD simulations, to characterize the morphology beyond the nano-scale. This theoretical approach enables us to comprehensively document the molecular-scale impact of the 1-CN molecules within the PM6:Y6 blends as a function of its concentration in solution, aspects that remain experimentally challenging to characterize. Specifically, regarding the “global” morphology, we examine the evolution in size, connectivity, and crystallinity of the PM6 and Y6 domains as a function of 1-CN concentration; regarding the local morphology, we describe the evolution of the PM6_PM6, Y6_Y6, and PM6_Y6 packing order, distance, and pattern as a function of 1-CN concentration; regarding the electronic properties, we discuss the evolution of the electron/hole transfer rates between adjacent Y6/PM6 molecules/chains, energetic distributions of the interfacial charge-transfer (CT) electronic states, and non-radiative recombination rates and related voltage losses from the CT states to the ground state, as well as the hole/electron transfer rates between PM6/Y6 and 1-CN moieties. Our results allow us to rationalize, at the molecular scale, how the 1-CN molecules systematically alter the photovoltaic parameters and thus the PCEs of PM6:Y6-based OSCs, and to point out the requirements in terms of the electronic properties of solvent additives to lead to further improvements in device efficiency.

Results and discussion

Evolution of the morphology at the global level

As a first step, we examine the evolutions in size, connectivity, and crystallinity of the Y6 and PM6 domains in the PM6:Y6:1-CN blends as a function of 1-CN concentration. It is worth pointing out that producing highly efficient OSCs usually calls for moderate domain sizes on the order of a few tens of nanometers and a large extent of domain connectivity and crystallinity.^9,15 Fig. 2 presents the morphologies of the PM6:Y6:1-CN blends at 1-CN concentrations in solution of 0%, 0.25%, 0.5%, and 1% v/v, respectively; these are obtained via solvent-evaporation simulations in the framework of CG-MD. Since there is a lack of experimental reports determining the amount of 1-CN molecules remaining in the active layer after solvent evaporation, here, we make the assumption that all are left within the PM6:Y6 blends. The critical aspect is that this procedure will not alter the trends obtained as a function of 1-CN concentration. In each blend, there exist four types of regions: Y6, PM6, mixed Y6:PM6, and 1-CN regions. As the concentration in 1-CN increases, the degree of 1-CN aggregation also increases while the volume fractions in Y6, PM6, and mixed regions decrease. The increase in the extent of 1-CN aggregation is seen in Fig. 2 to trim the Y6 as well as PM6 regions, which is expected to lead to smaller domain size and weaker domain connectivity. To confirm these points, we extracted the largest Y6 clusters in terms of number of Y6 molecules and PM6 clusters in terms of number of repeat units from these blends. Indeed, as the 1-CN concentration increases from 0%, 0.25%, 0.5%, to 1% v/v, the sizes of the largest Y6 and PM6 clusters both significantly reduce (see Table S1). This indicates that 1-CN can diminish the aggregation of Y6 molecules as well as of PM6 chains, which is consistent with transmission electron microscope (TEM) data showing that the use of 1-CN effectively suppresses the large and uneven phase separations leading to oversized domain sizes.^6,7 In this context, it is important to note that it has been widely reported that an excessive aggregation of Y6 molecules in the PM6:Y6 blend results in oversized Y6 domains and limits the OSC efficiency.^6,7,34–36

Fig. 2 In-plane slices of CG-MD simulated PM6:Y6:1-CN blends as well as their Y6 (blue), PM6 (red), mixed (yellow), and 1-CN (green) regions at 1-CN concentrations in solution of 0%, 0.25%, 0.5%, and 1% v/v, respectively. The percentages below each image represent the volume fractions of the various regions in the blends.

To characterize the intra-domain connectivity within Y6 and PM6 regions, we evaluated the ratios between the volumes spanned by the largest clusters (where a cluster is defined as a network of molecules/chains separated by at most 5 Å from their neighbors) and the overall volumes of the blends. The relevant spatial distributions of the largest Y6 and PM6 clusters in the CG-MD-simulated PM6:Y6:1-CN blends are provided in Fig. S1. As shown in Table S1 and Fig. S1, two main results are observed: (i) as the 1-CN concentration raises, the degree of connectivity in both Y6 and PM6 regions decreases; and (ii) at the 1-CN concentration of 1% v/v (corresponding in reality to a higher 1-CN concentration in solution due to our assumption that all 1-CN molecules remain within the PM6:Y6 blend), the connectivity within the Y6 region is significantly disrupted, which is expected to lead to a decrease in carrier mobility and an increase in charge recombination and thus a lower PCE. We note that changing the threshold (here, ∼5 Å) for the intermolecular/interchain distance used to extract the Y6/PM6 clusters would alter the absolute values but not the trends illustrated in Table S1.

In order to examine the evolution of crystallinity within the Y6, PM6, and mixed regions as a function of 1-CN concentration, the radial distribution functions (RDFs), g(r), were analyzed for the Y6_Y6 backbones (i.e., RDFs between Y6 backbones), PM6_PM6 backbones (i.e., RDFs between PM6 backbones), and PM6_Y6 backbones (i.e., RDFs between PM6 and Y6 backbones) in the CG-MD-simulated PM6:Y6:1-CN blends. We recall that the RDF measures how the relative density of particles (such as atoms or molecules) varies as a function of the distance away from a reference particle; if ρ_global is the average number density of particles, then the local averaged density at a distance r is ρ_local = ρ_globalg(r).^37–39 The first peak represents the relative density of nearest-neighbor stacking, the second peak indicates that of next nearest-neighbor stacking, and so forth. Therefore, a larger g(r) peak value points to a higher extent of crystallinity. It is clear from Fig. 3a, b, and c that, as the 1-CN concentration increases, the g(r) peak values increase in all the Y6_Y6, PM6_PM6, and PM6_Y6 cases; this suggests an enhancement in the extent of crystallinity or stacking order in Y6, PM6, and mixed regions upon addition of 1-CN molecules. We note that it is inappropriate in this instance to use the relative heights of the g(r) peaks to identify which blend has the larger Y6_Y6, PM6_PM6, or PM6_Y6 stacking density since the ρ_global value differs among PM6:Y6:1-CN blends with different 1-CN concentrations.

Fig. 3 Radial distribution functions for (a and d) Y6_Y6 backbones, (b and e) PM6_PM6 backbones, and (c and f) PM6_Y6 backbones in the CG-MD-simulated/AA-MD-simulated PM6:Y6:1-CN blends.

Evolution of the morphology at the local (nanoscale) level

We now turn to the evolution of the local morphology by examining the Y6_Y6, PM6_PM6, and PM6_Y6 packing order, distance, and pattern for the AA-MD-simulated PM6:Y6:1-CN blends with different 1-CN concentrations. Considering that our CG-MD simulations have shown that an excess in 1-CN molecules segments the Y6 and PM6 regions into very small and disconnected domains, which is expected to break the charge-transport pathways, we only focus here on 1-CN concentrations of 0%, 0.25%, and 0.5%. Fig. 3d, e, and f display the RDFs for the Y6_Y6, PM6_PM6, and PM6_Y6 backbones in these blends. What can be observed is that adding 1-CN molecules (i.e., from 0% to 0.25% and 0.5% v/v) increases the number of g(r) peaks for all the Y6_Y6 (i.e., from 2 to 3), PM6_PM6 (i.e., from 1 to 2), and PM6_Y6 (i.e., from 1 to 2) cases. In other words, the 1-CN molecules contribute to enhance the Y6_Y6, PM6_PM6, and PM6_Y6 packing order, as also confirmed when examining the CG-MD-simulated PM6:Y6:1-CN blends (see Fig. 3a, b, and c). In addition, along with an increase in 1-CN concentration from 0% to 0.25% and then 0.5% v/v, the Y6_Y6, PM6_PM6, and PM6_Y6 backbone–backbone packing distances decrease from 4.25 to 4.15 and then 4.05 Å, from 4.95 to 4.65 and then 4.35 Å, and from 4.65 to 4.35 and then 4.25 Å, respectively, as supported by comparing the first g(r) peak positions displayed in Fig. S2.

To understand these findings, we extracted all the Y6_Y6, PM6_PM6, and PM6_Y6 π–π packing pairs from the AA-MD-simulated PM6:Y6:1-CN blends and computed their averaged respective interaction energies. Here, a π–π packing pair is defined as consisting of two backbones having their atoms directly interacting face-to-face within a distance of 5 Å of each other. As summarized in Table S2, in going from 0% to 0.25% and then 0.5% v/v, the interactions among the Y6_Y6, PM6_PM6, and PM6_Y6 π–π packing pairs all become stronger (i.e., from −16.1 ± 6.2 to −20.1 ± 6.8 and then −20.3 ± 6.8 kcal mol⁻¹, from −25.9 ± 7.4 to −27.4 ± 10.8 and then −28.3 ± 11.4 kcal mol⁻¹, and from −21.9 ± 8.7 to −25.2 ± 10.1 and then −25.9 ± 10.2 kcal mol⁻¹, respectively). The increase in the Y6_Y6, PM6_PM6, and PM6_Y6 interaction strength goes hand in hand with the increase [decrease] in the Y6_Y6, PM6_PM6, and PM6_Y6 packing order [distance].

Considering that the Y6 backbone is composed of core and end moieties and the PM6 backbone, of D (electron-rich) and A (electron-poor) moieties (see Fig. S3), there typically exist three Y6_Y6 [i.e., end_end (end moieties on top of one another), end_core (end moiety on top of core moiety), and core_core (core moieties on top of one another)] π–π packing configurations, three PM6_PM6 [i.e., A_A (A moieties on top of one another), A_D (A moiety on top of D moiety), and D_D (D moieties on top of one another)] π–π packing configurations, and four PM6_Y6 [i.e., A_end (end moiety on top of A moiety), D_end (end moiety on top of D moiety), D_core (core moiety on top of D moiety), and A_core (core moiety on top of A moiety)] π–π packing configurations. To explore the evolution of the Y6_Y6, PM6_PM6, and PM6_Y6 packing patterns as a function of the 1-CN concentration, Fig. 4 displays the partial RDFs for the Y6_Y6, PM6_PM6, and PM6_Y6 backbones in the AA-MD-simulated PM6:Y6:1-CN blends with different 1-CN concentrations. Interestingly, the 1-CN molecules have negligible effect on the preferential Y6_Y6, PM6_PM6, and PM6_Y6 π–π packing orders [i.e., end_end > end_core > core_core; A_A > A_D > D_D; and A_end > D_end > D_core > A_core, respectively] (see Fig. 4), as also confirmed by the partial RDFs for the Y6_Y6, PM6_PM6, and PM6_Y6 backbones in the CG-MD-simulated PM6:Y6:1-CN blends with different 1-CN concentrations. These packing orders result mainly from: (i) the steric hindrance induced by the alkyl side chains on the Y6 core moieties (see Fig. S3a); (ii) the torsion of the fluorothiophenes (in the 2-(2-ethylhexyl)fluorothiophene side chains) with respect to the D moieties of the PM6 backbones (in addition, while the PM6 D moieties carry branched alkyl side chains on both sides, the PM6 A moieties carry them on only one side, which leaves the other side of the PM6 A moieties more open for interaction with Y6, see Fig. S3b); and (iii) the larger sizes of the Y6 core moieties and PM6 D moieties (see Fig. S3), which bring about stronger D_core than A_core interactions and thus leads to a preferential packing order of D_core > A_core.

Fig. 4 Partial radial distribution functions for the Y6_Y6, PM6_PM6, and PM6_Y6 backbones in the AA-MD-simulated PM6:Y6:1-CN blends.

Evolution of the electronic properties related to the nanoscale morphology

Having in hands all the Y6_Y6, PM6_PM6, and PM6_Y6 π–π packing pairs extracted from the AA-MD-simulated PM6:Y6:1-CN blends, we now examine the evolution of the electronic properties as a function of 1-CN concentration. Table S3 collects the intermolecular electron transfer rates, k_e, between two neighboring Y6 molecules (evaluated on the basis of the semi-classical Marcus theory)⁴⁰ and their related parameters for the Y6_Y6 π–π packing pairs. As the 1-CN concentration increases from 0% v/v to 0.25% and 0.5% v/v, the k_e values slightly increase from 2.0 × 10¹² s⁻¹ to 3.4 × 10¹² and 3.7 × 10¹² s⁻¹ (for transfer from a Y6 molecule with a smaller electron affinity, EA, to one with a larger EA) and from 2.8 × 10¹¹ s⁻¹ to 5.0 × 10¹¹ and 4.9 × 10¹¹ s⁻¹ (for transfer from a Y6 with a larger EA to one with a smaller EA); this is mainly attributable to the increase in electronic couplings between the initial and final states (comparing the averaged parameters in Table S3). To further confirm these results, Table S4 summarizes the proportions of Y6_Y6 π–π packing pairs corresponding to various orders of magnitude in these rates. At each of the higher levels of electron transfer rates (i.e., ≥10¹³, ≥10¹², ≥10¹¹, ≥10¹⁰, and ≥10⁹), there basically exists a larger proportion of pairs in the case of the PM6:Y6:1-CN blends with 1-CN concentrations of 0.25% and 0.5% v/v.

Table S5 lists the interchain hole transfer rates, k_h, and the related microscopic parameters for the PM6_PM6 π–π packing pairs extracted from these AA-MD-simulations. Similar to what we found regarding the electron transfer rates among Y6 molecules, the k_h values increase slightly from 1.4 × 10¹² to 1.8 × 10¹² and then to 1.9 × 10¹² s⁻¹ (for transfer from a PM6 segment with a smaller ionization potential, IP, to one with a larger IP) and from 7.4 × 10¹⁰ to 8.0 × 10¹⁰ and then 8.4 × 10¹⁰ s⁻¹ (for transfer from a PM6 segment with a larger IP to one with a smaller IP), as a function of an increase in 1-CN concentration from 0% to 0.25% and then 0.5% v/v. To further validate these results, Table S6 collects the proportions of PM6_PM6 π–π packing pairs corresponding to various orders of magnitude in these rates. At each of the higher levels of hole transfer rates (i.e., ≥10¹³, ≥10¹², ≥10¹¹, ≥10¹⁰, and ≥10⁹), the proportion of pairs follows the order 0.5% v/v > 0.25% v/v > 0% v/v.

We now discuss the correlation of the interfacial PM6_Y6 π–π packing patterns with the characteristics of the lowest singlet CT electronic states, which are the key intermediate states in the charge generation process.⁴¹ Fig. 5a presents the energetic distributions of the CT states for the AA-MD-simulated PM6:Y6:1-CN blends with different 1-CN concentrations. In going from 0% to 0.25% and then to 0.5% v/v, there occurs a decrease in both the average energies of the CT states, E^avg_CT, and their standard deviations, σ: from 1.636 to 1.621 and then 1.619 eV for the former and from 0.124 to 0.114 and then 0.105 eV for the latter. We note that increasing the number of extracted the donor–acceptor packing pairs and modifying the DFT basis set (for instance, comparing 6-31G** and 6-31G*) do not change the trends in the evolution of both E^avg_CT and σ. The decrease in σ can be ascribed to the increase in the PM6_Y6 π–π packing order (as discussed above), which translates into lower voltage loss related to interfacial disorder.^42–44 It has been reported that the CT-state energy (E_CT) can be evaluated, to a first approximation, as the sum of the transport gap energy (E_gap) and the hole–electron electrostatic interaction energy (E_el, for which a negative value means attraction):⁴⁴

E_CT = E_gap + E_el = IP_donor − EA_acceptor + E_el,

where IP_donor and EA_acceptor are the ionization potential of the donor (here, PM6) and the electron affinity of the acceptor (here, Y6), respectively. Fig. 5b shows the evolution of E_CT, E_gap, and E_el as a function of 1-CN concentration. Clearly, E_CT in this instance has a stronger correlation with E_el than E_gap (i.e., E_CT decreases as E_el becomes increasingly negative). In other words, as the 1-CN concentration increases, the PM6_Y6 packing distances decrease (see Fig. S2c), which leads to a stronger hole–electron electrostatic attraction and thus a lower E_CT value.

Fig. 5 AA-MD-simulated PM6:Y6:1-CN blends with 1-CN concentration of 0%, 0.25%, and 0.5% v/v: (a) normalized energetic distributions of the CT states, fitted with Gaussian functions; (b) CT-state energies E_CT, transport gap energies E_gap, and hole–electron electrostatic interaction energies E_el; (c) non-radiative recombination rates k_nr; and (d) non-radiative voltage losses ΔV_nr as a function of the considered quantum-mechanical component of the reorganization energy λ_qm.

A recent work by Pratik et al. has demonstrated that there also exist PM6_Y6 π–π packing pairs having a second Y6 molecule positioned near, and interacting significantly with the initial Y6 molecule;⁴⁵ the presence of the second Y6 molecule is expected to decrease both the E_el value (as the hole wavefunction can become more delocalized) and the E_gap value (as the Y6_Y6 interaction will reduce the local IP). To estimate more quantitatively the impact that the second Y6 molecule can have, we randomly extracted three PM6_Y6_Y6 π–π packing complexes from the AA-MD-simulated PM6:Y6 blend and compared their averaged E_CT, E_gap, and E_el values with those of the corresponding PM6_Y6 pairs. The results show that the presence of the second Y6 molecule decreases the hole–electron electrostatic attraction strength (i.e., a less negative E_el value) to a lesser extent than the transport gap, which overall leads to a slightly lower E_CT value (see Table S7).

Fig. 5c displays in the case of these AA-MD-simulated blends the non-radiative recombination rates from the lowest singlet CT state to the ground state, k_nr, as a function of the quantum-mechanical component of the reorganization energy, λ_qm. It has been widely demonstrated that in the instance of high CT-state energies, the quantum-mechanical component (due to coupling to high-frequency vibration modes) of the reorganization energy becomes especially important.^46,47 Thus, we turned to the Marcus–Levich–Jortner model to consider the non-radiative recombination rates.^46,48,49 As shown in Fig. 5c, the blends in the presence of 1-CN (0.25% and 0.5% v/v) have consistently higher k_nr values across the whole range of λ_qm considered, which can be assigned to their lower average CT-state energies and higher average electronic couplings by comparing the parameters listed in Table S8. The implication of this finding is that higher non-radiative voltage losses, ΔV_nr, are expected in the AA-MD-simulated PM6:Y6:1-CN blends with 1-CN concentration of 0.25% and 0.5% v/v. To confirm this point, the ΔV_nr values were estimated by considering:^50,51

where EQE_EL represents the external quantum efficiency of electroluminescence; k_B, the Boltzmann constant; T, the temperature (298.15 K); q, the electron charge; p_e, the probability for a photon generated by radiative recombination to escape from the device; and k_r, the radiative recombination rates from the lowest singlet CT state to the ground state (the latter rates are evaluated by the Einstein relation,⁵² with the relevant parameters provided in Table S9). Fig. 5d displays the ΔV_nr values as a function of λ_qm for the PM6:Y6:1-CN blends with different 1-CN concentrations; here, the p_e value is taken to be 0.25, which we note lies in the range of the largest values of ∼0.2–0.3 measured for the light-outcoupling coefficient in organic light-emitting devices.^38,53 Over the whole range of λ_qm values we considered, the ΔV_nr values are indeed higher for the PM6:Y6:1-CN blends with 1-CN concentration of 0.25% and 0.5% v/v; this is in accord with the experimental findings that the PM6:Y6-based OSCs fabricated by adding 0.5% v/v 1-CN in chloroform solution possess slightly higher ΔV_nr values (250 meV vs. 240 meV)¹⁹ and lower V_OC values (0.82–0.84 V vs. 0.85 V) with respect to the solar cells fabricated without additive.^4,5,11,19

Key requirements regarding the electronic properties of residual solvent additives

Finally, we highlight the key requirements in terms of the electronic properties of residual solvent additives. As an initial step, we identify to what extent the 1-CN molecules interact with PM6 chains and Y6 molecules. To do so, we extracted the PM6_1-CN π–π packing pairs from the AA-MD-simulated PM6:Y6:1-CN blend with the 1-CN concentration of 0.5% v/v that is widely adopted in experimental studies.^{4–7,11,19,23,28} We obtain 207 PM6_1-CN π–π packing pairs, while following the same procedure, we obtain 375 Y6_1-CN π–π packing pairs. Thus, the 1-CN molecules are present in both PM6 and Y6 domains, but more so in the latter. Our analyses of the RDFs for the PM6/Y6_1-CN backbones (i.e., RDFs between PM6/Y6 backbones and 1-CN backbones) in the AA-MD-simulated PM6:Y6:1-CN blend (0.5% v/v) further confirm this point by showing a clear first g(r) peak around ∼5 Å for both the PM6_1-CN and Y6_1-CN cases but a relatively higher first g(r) peak for the latter (see Fig. S5).

To place this result on thermodynamical grounds, we evaluated the Flory–Huggins interaction parameters, χ, of the PM6:1-CN and Y6:1-CN blends, using the equation:^54,55

where K is a positive constant and γ₁ and γ₂ are the surface energies of components 1 and 2, respectively. We emphasize that the χ parameter offers an overall measure of the intermolecular interactions present in a blend; a smaller χ value indicates a larger extent of mixing between two components. Based on the experimental surface energies of PM6 (19.12 mN m⁻¹), Y6 (19.88 mN m⁻¹), and 1-CN (20.65 mN m⁻¹),⁸ the χ value of the Y6:1-CN blend is estimated to be smaller than that of the PM6:1-CN blend (i.e., 0.561 K vs. 0.602 K), which is consistent with the greater number of 1-CN molecules interacting with Y6 in the AA-MD-simulated PM6:Y6:1-CN blend. It is worth noting that there exist large deviations among the surface energies coming from different experimental measurements;^8,30,56 to ensure the reliability of the χ calculations, we adopted the values of the PM6, Y6, and 1-CN surface energies coming from the same work.⁸

To address the question of whether the 1-CN molecules participate in the charge transport processes taking place in the PM6:Y6:1-CN blend, we examined the charge transfer rates between PM6/Y6 and 1-CN, which involve the transfer processes from PM6/Y6 to 1-CN as well as their reverse. In the framework of the semi-classical Marcus theory,⁴⁰ the hole transfer rates, k_h, for the 207 PM6_1-CN π–π packing pairs are calculated to be in the range of 8 × 10⁷ s⁻¹ to 4 × 10¹⁴ s⁻¹ (from PM6 to 1-CN) and 10⁻⁴ s⁻¹ to 10⁸ s⁻¹ (from 1-CN to PM6), respectively, see Table S10; the electron transfer rates, k_e, for the 375 Y6_1-CN π–π packing pairs are essentially 0 for transfer from Y6 to 1-CN and in the range of 10⁻³⁸ s⁻¹ to 4 × 10⁻⁸ s⁻¹ for transfer from 1-CN to Y6, see Table S11. Considering that the k_h values (from 1-CN to PM6) and k_e values (from Y6 to 1-CN and vice versa) are much lower than those between PM6 chains and Y6 molecules (as discussed above), it can be concluded that the 1-CN molecules can hardly get involved in the charge transport processes occurring in the PM6:Y6:1-CN blend. As shown in Tables S10 and S11, this result can be attributed to the large energy difference in IPs between 1-CN and PM6 (i.e., ∼0.52 eV) and in EAs between Y6 and 1-CN (i.e., ∼1.95 eV).

The above finding motivated us to examine: (i) what would be the requirements that the 1-CN IP or EA would need to fulfill in order for the additive to get involved in the charge transport processes in the PM6:Y6:1-CN blend (i.e., leading to k_h [k_e] values at least comparable to those between PM6 chains [Y6 molecules]); and (ii) whether 1-CN could take part in the charge transport processes in blends based on derivatives of PM6:Y6. Fig. 6 answers question (i) by showing the evaluated k_h [k_e] values between PM6 [Y6] and 1-CN as a function of modulating the IP_PM6–IP_1-CN [EA_Y6–EA_1-CN] values; here, the calculations of k_h and k_e are based on the electronic couplings and reorganization energies averaged over the corresponding data in Tables S10 and S11, respectively. It is clear from Fig. 6 that, for the k_h [k_e] values to be comparable to those between PM6 chains [Y6 molecules], the |IP_PM6–IP_1-CN| [|EA_Y6–EA_1-CN|] difference would need to be less than ∼0.15 eV [∼0.10 eV]; note that k_h [k_e] values between PM6 chains [Y6 molecules] are on the order of 10¹⁰ s⁻¹ [10¹¹ s⁻¹]. Since the actual |IP_PM6–IP_1-CN| and |EA_Y6–EA_1-CN| values are as large as ∼0.52 eV and ∼1.95 eV, respectively, the 1-CN molecules, as concluded above, cannot participate effectively in the charge transport processes in the PM6:Y6:1-CN blend.

Fig. 6 (a) Evaluated hole transfer rates, k_h, between PM6 and 1-CN, as a function of the modulation of the IP_PM6–IP_1-CN value; (b) electron transfer rates, k_e, between Y6 and 1-CN, as a function of the modulation of the EA_Y6–EA_1-CN value.

With the requirements derived from Fig. 6 in mind, we turn to answer the question (ii). Fig. S6a displays the chemical structures of widely used PM6 and Y6 derivatives, while Fig. S6b provides their IP and EA values with respect to those of 1-CN. It is clear that the |IP_donor–EA_1-CN| and |EA_acceptor–EA_1-CN| values for the relevant donors and acceptors are much greater than ∼0.15 eV and ∼0.10 eV, respectively, which suggests that the 1-CN residues would also hardly play any significant role in the charge transport processes in the blends based on PM6:Y6 derivatives (at least for those shown in Fig. S6a).

Overall, our results highlight that the 1-CN molecules can have intermolecular interactions with both PM6 chains and Y6 molecules, but do not partake in the charge transport processes in the PM6:Y6:1-CN blends. Therefore, in order to further improve the PCEs of OSCs based on the PM6:Y6 blend or its derivatives, this calls for the design of novel solvent additives to not only optimize the active-layer morphology but also finely tune the charge transport processes therein (e.g., by creating new charge transport pathways to simultaneously enhance and balance carrier mobilities).

Conclusions

Examining the morphological and electronic properties of PM6:Y6:1-CN blends as a function of the additive concentration in solution, through a combination of long-range corrected DFT calculations, AA-MD simulations, and CG-MD simulations, has enabled us to draw a comprehensive picture of how the 1-CN molecules impact the PM6:Y6 blends and their photovoltaic parameters at the molecular scale. The results allowed us to point out some specific requirements regarding the electronic properties of solvent additives to further improve device efficiency. The main conclusions of our work are as follows:

(i) As the 1-CN concentration increases, the crystallinity of the Y6, PM6, and their mixed domains improves; in addition, the extent of 1-CN aggregation increases, which decreases the size of the Y6 and PM6 domains and their inner connectivity.

(ii) As the 1-CN concentration increases, the interaction strengths among different components increase, which decreases the Y6_Y6, PM6_PM6, and PM6_Y6 π–π packing distances and increases their packing order.

(iii) The decrease in the Y6_Y6 and PM6_PM6 π–π packing distances leads to an increase in the electronic couplings relevant for the charge transfer processes, and thus improves electron/hole transfer rates between adjacent Y6 molecules/PM6 chains. Furthermore, the decrease in the PM6_Y6 π–π packing distances increases the hole–electron electrostatic attraction as well as the electronic coupling relevant for non-radiative recombination process; the former decreases the CT-state energies, which together with the latter bring about higher non-radiative recombination rates and voltage losses, which is consistent with the experimental data.

(iv) The 1-CN molecules can interact with both PM6 chains and Y6 molecules and be present in both domains but can hardly participate in the charge transport processes taking place in the PM6:Y6 blend. Similar results are expected for blends based on PM6:Y6 derivatives. To get involved in the hole [electron] transfer processes, the IP [EA] of the residue molecules should be as close as possible (within 0.15 eV [0.10 eV]) to that of PM6 [Y6].

Based on our results, the factors that most account for the improved photovoltaic parameters and thus the PCEs of PM6:Y6:1-CN-based OSCs at a specific 1-CN concentration can be summarized as follows: (i) preventing excessive aggregation of Y6 molecules/PM6 chains (i.e., providing more appropriate domain sizes), which contributes to an increase in J_SC by producing a greater PM6:Y6 interfacial area and more efficient exciton dissociation and by shortening the distance for exciton diffusion to the PM6:Y6 interfaces; (ii) improving the crystallinity within the Y6, PM6, and their mixed domains, which effectively increases carrier mobility and thus J_SC and FF; and (iii) decreasing the Y6_Y6 and PM6_PM6 packing distances, which improves carrier mobility, J_SC, and FF by increasing the electronic couplings related to the charge transfer processes and boosting both electron and hole transfer rates. These features overall lead to higher carrier mobility, short-circuit current density, and fill factor. To further improve the OSC efficiency, efforts are needed to design novel solvent additives that could effectively participate in the charge transport processes, enhance the size of molecular networks within the domains, and increase the interfacial CT-state energies, while keeping the above advantages.

Computational methodology

All-atom molecular dynamics simulations

The AA-MD simulations were carried out with the LAMMPS package⁵⁷ and the optimized potentials for liquid simulations-all atom (OPLS-AA) force field.^37,38,58,59 To accurately describe the intramolecular and intermolecular interactions for the PM6, Y6, and 1-CN components, the OPLS-AA force field was reparametrized based on the results of long-range corrected DFT calculations: (i) the atomic partial charges were derived by fitting the electrostatic potential (ESP) from DFT calculations carried out at the ωB97XD/cc-PVTZ level of theory; (ii) the bond lengths and angles were taken from the geometries optimized at the ωB97XD/6-31G(d,p) level of theory, with the harmonic force constants kept as they are; and (iii) the parameters for the inter-ring dihedrals along the PM6 chains or Y6 molecules were fitted on the basis of their torsion potentials evaluated at the ωB97XD/6-31G(d,p) level of theory. Here, the relevant DFT calculations were performed with the Gaussian 16 C.01 package⁶⁰ and the ω range-separation parameters were optimized in the gas phase, i.e., for isolated molecules.

The initial models for the PM6:Y6:1-CN blends with 1-CN concentration of 0%, 0.25%, 0.5% v/v in solution were built by randomly placing the PM6 chains (each chain consisting of 8 repeat units), Y6 molecules, and 1-CN molecules in three cubic boxes with a low density of 0.02 g cm⁻³. The boxes that correspond to the cases of 0% v/v, 0.25% v/v, 0.5% v/v include 30 PM6 chains and 242 Y6 molecules together with 0, 720, and 1440 1-CN molecules, respectively. The AA-MD simulations were first carried out with the NPT (constant number of molecules, pressure, and temperature) ensemble for 30 ns at a temperature of 500 K and a pressure of 1 atm. Then, the three blends were rapidly cooled down from 500 K to 298.15 K, with a cooling rate of ∼50 K ns⁻¹. To simulate the experimental annealing process, the AA-MD simulations were finally performed (i) at 383.15 K for 20 ns; (ii) from 383.15 K to 298.15 K with a slow cooling rate of ∼2 K ns⁻¹; and (iii) at 298.15 K for 10 ns; we note that 383.15 K is the experimental annealing temperature.^4,16,28 A cutoff of 12 Å was used for the summation of the van der Waals interactions and the particle–particle particle–mesh (PPPM) solver for the long-range Coulomb interactions. The Verlet integrator was considered with a timestep of 1 fs and the Nosé–Hoover thermostat/barostat was employed for temperature/pressure control.

It is worth noting that varying the numbers of PM6 chains, Y6 molecules, and 1-CN molecules proportionally or altering the annealing temperature within a reasonable range (for instance, from an experimental temperature of 383.15 K (ref. 4, 16 and 28) to another experimental temperature of 373.15 K (ref. 10, 19 and 31)) might change the absolute values of the data but not any of the conclusions of our work.

Coarse-grained molecular dynamics simulations

The CG-MD simulations were performed using the GROMACS package and the Martini 3 CG force field.⁶¹ To accurately capture the intra- and intermolecular interactions for the PM6, Y6, and 1-CN components, the Martini 3 force field was refined based on the results from our AA-MD simulations: (i) the selection of CG bead types for different segments was based on the small-molecule building-block table provided by the Martini 3 force field,⁶² which specifies the CG bead types and the related chemical structural units. For chemical structural units that are not listed in that table, the optimal bead types were selected by comparing the calculated transfer free energies with the experimental measurements and by comparing potentials of mean force (PMF) from the CG-MD simulations to those from the AA-MD simulations (the CG representations of PM6, Y6, and 1-CN are shown in Fig. S7); (ii) the bond lengths and angles for beads were derived from the AA-MD simulations. Specifically, the trajectories were first generated via 50 ns AA-MD simulations and then mapped to the CG bead trajectories according to the established bead mapping scheme.⁶³ The distributions of bond lengths and angles obtained in these mapped trajectories were subsequently used to determine the equilibrium bond lengths, angles, and harmonic force constants for the Martini 3 force field; and (iii) the dihedral parameters between cyclic structures formed by the CG beads were again derived from the mapped CG trajectories. The distributions of bead dihedrals obtained from these trajectories were fitted by using the Ryckaert–Bellemans functional form ref. 63 to parameterize proper dihedrals within the Martini 3 force field. This procedure ensured consistency between dihedral distributions from the CG-MD simulations and those from the mapped AA-MD trajectories.

Four systems were built by adding 1-CN molecules to PM6:Y6 solutions at concentrations of 0%, 0.25%, 0.5%, and 1% v/v, respectively, with chloroform taken as solvent, the weight ratio between PM6 and Y6 set at well 1 : 1.2, and the PM6:Y6 concentration considered as ∼20 mg ml⁻¹. Here, a PM6 chain consists of 30 repeat units, which corresponds to a number-average molecular weight of ∼36.66 kDa. The systems corresponding to 0%, 0.25%, 0.5% and 1% v/v 1-CN include 92 PM6 chains, 2790 Y6 molecules, and 4 574 270 CF molecules, together with 0, 11 088, 22 176, and 44 352 1-CN molecules, respectively. The solvent-evaporation simulations for each system were performed following the protocol proposed by Alessandri and co-workers.^64,65 During the simulations, the temperature [pressure] was controlled using the velocity-rescaling thermostat [Parrinello–Rahman barostat]. A cutoff of 12 Å was applied for the van der Waals interactions and the Coulomb interactions were treated using the reaction-field method. Given that the 1-CN molecules have been demonstrated to retard the solvent-evaporation process, the time interval between two solvent-removal steps was set to 2 ns for system without 1-CN and to 4 ns for systems with 1-CN; at each removal step, 1.25% of the chloroform molecules were removed. The integration time step was set to 20 fs.

The final dimensions for the four CG-MD simulated PM6:Y6:1-CN blends are as follows: for the 0% case, 31.2 × 31.2 × 8.5 nm³; for the 0.25% case, 31.8 × 31.8 × 10.2 nm³; for the 0.5% case, 29.8 × 29.8 × 13.9 nm³; and for the 1% case, 31.2 × 31.2 × 17.1 nm³.

Density functional theory calculations

The intermolecular backbone–backbone interaction energies for the Y6_Y6, PM6_PM6, and PM6_Y6 π–π packing pairs were calculated at the ωB97XD/6-31G(d,p) level of theory, where all pair geometries were kept as those extracted from the AA-MD-simulated blends; all side chains were replaced with methyl groups. To avoid the overstabilization of these energies induced by basis set superposition error, the counterpoise correction method proposed by Boys and Bernardi was employed.⁶⁶

The semi-classical Marcus theory was used to evaluate the charge transfer rate, k_i (i = e or h for electron or hole), between two molecules in a π–π packing pair:⁴⁰

where ħ represents the reduced Planck constant; V_i, the electronic coupling between the initial and final states; λ_i, the reorganization energy; and ΔE_i, the energy difference between the initial and final states. V_i between the initial and final states was approximated as the coupling between the frontier molecular orbitals (i.e., lowest unoccupied molecular orbital for electron and highest occupied molecular orbital for hole) of the two molecules in the π–π packing pair; V_i was estimated by the fragment orbital method.⁶⁷ ΔE_i between the initial and final states was estimated as the energy difference between ionization potential (IP) or electron affinity (EA) (i.e., EA for electron and IP for hole) of the two molecules in the π–π packing pair. The λ_i reorganization energy consists of intra- (λ_intra) and intermolecular (λ_inter) contributions. The former was evaluated from the adiabatic potential energy surfaces of the molecular states involved in the charge transfer process; the latter was set to 0.1 eV (a reasonable value for extended π-conjugated systems),^68,69 considering that there are currently no straightforward models to accurately estimate the λ_inter value in solid-state environments. The DFT calculations were performed at the PCM-tuned-ωB97XD/6-31G(d,p) level of theory (here, “PCM-tuned” means that the ω range-separation parameter was optimized within the polarizable continuum model and considering a typical dielectric constant of 3.5).^43,70

The Marcus–Levich–Jortner model was adopted to examine the non-radiative recombination rates, k_nr, from the CT states to the ground state:^46,48,49

where V_el denotes the electronic coupling between the CT state and the ground state; S_qm, the Huang–Rhys parameter with S_qm = λ_qm/ħω_qm; λ_qm, the quantum-mechanical contribution to the reorganization energy, described by means of an effective vibrational mode ħω_qm; λ_c, the classical contribution to the reorganization energy; and E^ad_CT, the adiabatic (relaxed) CT energy. V_el between the CT state and the ground state was estimated in the framework of the generalized Mulliken–Hush approach.⁷¹ E^ad_CT was evaluated as: E^ad_CT = E^avg_CT − λ, where E^avg_CT is the averaged vertical transition energies of the CT states computed (via time-dependent density functional theory, TD-DFT) for PM6_Y6 π–π packing pairs at their ground-state geometries as extracted from the AA-MD-simulated blends. The λ_intra value was estimated from the adiabatic potential energy surfaces of the molecular states involved in the nonradiative recombination process; λ_inter was again assumed to be 0.1 eV.^68,69 Considering that the PM6_Y6 π–π packing pairs extracted from the AA-MD-simulated blends are too large for full geometry optimizations and frequency calculations of the excited states, we assumed a typical ω_qm value of 1200 cm⁻¹.^38,72 We emphasize that altering this value within a reasonable range will not change any of the conclusions. Moreover, considering that the exact partition of the overall reorganization energy into classical and quantum contributions is unknown, we considered the evolution of the k_nr values as a function of λ_qm, where λ_c was taken as λ_c = λ − λ_qm. Here, the DFT and TD-DFT calculations were carried out at the PCM-tuned-ωB97XD/6-31G(d,p) level of theory.

The Einstein coefficient relation was used to evaluate the radiative recombination rates, k_r, from the CT states to the ground state:⁵¹

where ε₀ represents the vacuum permittivity; c, the speed of light; and μ_CT→s0, the transition dipole moment between the CT state and the ground state.

Author contributions

Ziwen Yu: writing – original draft, validation, methodology, formal analysis. Tongrui Zhang: writing – original draft, validation, methodology, formal analysis. Yizhou Li: methodology. Gao-Feng Han: writing – review. Xing-You Lang: writing – review. Jean-Luc Brédas: writing – review & editing, funding acquisition. Tonghui Wang: writing – original draft, writing – review & editing, methodology, formal analysis, funding acquisition. Qing Jiang: writing – review & editing, funding acquisition.

Conflicts of interest

The authors declare no competing interests.

Data availability

The data supporting this article have been included as part of the SI. Supplementary information is available. See DOI: https://doi.org/10.1039/d5el00117j.

Acknowledgements

The work at Jilin University was financially supported by the National Key R&D Program of China (No. 2023YFB3003001), National Natural Science Foundation of China (No. 52130101, 52271217), “Xiaomi Young Scholar” Project, and the fund of “World-class Universities and World-class Disciplines”, Ministry of Education, China. The work at the University of Arizona was funded by the Office of Naval Research, Award No. N00014-24-1-2114.

Notes and references

1. N. S. Sariciftci, L. Smilowitz, A. J. Heeger and F. Wudl, Science, 1992, 258, 1474–1476.
2. J. J. M. Halls, C. A. Walsh, N. C. Greenham, E. A. Marseglia, R. H. Friend, S. C. Moratti and A. B. Holmes, Nature, 1995, 376, 498–500.
3. G. Yu, J. Gao, J. C. Hummelen, F. Wudl and A. J. Heeger, Science, 1995, 270, 1789–1791.
4. J. Yuan, Y. Zhang, L. Zhou, G. Zhang, H.-L. Yip, T.-K. Lau, X. Lu, C. Zhu, H. Peng, P. A. Johnson, M. Leclerc, Y. Cao, J. Ulanski, Y. Li and Y. Zou, Joule, 2019, 3, 1140–1151.
5. J. Lv, H. Tang, J. Huang, C. Yan, K. Liu, Q. Yang, D. Hu, R. Singh, J. Lee, S. Lu, G. Li and Z. Kan, Energy Environ. Sci., 2021, 14, 3044–3052.
6. H. Fan, H. Yang, Y. Wu, O. Yildiz, X. Zhu, T. Marszalek, P. W. M. Blom, C. Cui and Y. Li, Adv. Funct. Mater., 2021, 31, 2103944.
7. C. Fan, H. Yang, Q. Zhang, S. Bao, H. Fan, X. Zhu, C. Cui and Y. Li, Sci. China: Chem., 2021, 64, 2017–2024.
8. X. He, C. C. S. Chan, J. Kim, H. Liu, C. Su, U. Jeng, H. Su, X. Lu, K. S. Wong and W. C. H. Choy, Small Methods, 2022, 6, 2101475.
9. J. Yi, G. Zhang, H. Yu and H. Yan, Nat. Rev. Mater., 2023, 9, 46–62.
10. L. Chen, J. Yi, R. Ma, L. Ding, T. A. Dela Peña, H. Liu, J. Chen, C. Zhang, C. Zhao, W. Lu, Q. Wei, B. Zhao, H. Hu, J. Wu, Z. Ma, X. Lu, M. Li, G. Zhang, G. Li and H. Yan, Adv. Mater., 2023, 35, 2301231.
11. X. Yang, B. Li, X. Zhang, S. Li, Q. Zhang, L. Yuan, D. Ko, W. Ma and J. Yuan, Adv. Mater., 2023, 35, 2301604.
12. F. Sun, X. Zheng, T. Hu, J. Wu, M. Wan, Y. Xiao, T. Cong, Y. Li, B. Xiao, J. Shan, E. Wang, X. Wang and R. Yang, Energy Environ. Sci., 2024, 17, 1916–1930.
13. L. Zhu, M. Zhang, G. Zhou, Z. Wang, W. Zhong, J. Zhuang, Z. Zhou, X. Gao, L. Kan, B. Hao, F. Han, R. Zeng, X. Xue, S. Xu, H. Jing, B. Xiao, H. Zhu, Y. Zhang and F. Liu, Joule, 2024, 8, 3153–3168.
14. L. Wang, C. Chen, Y. Fu, C. Guo, D. Li, J. Cheng, W. Sun, Z. Gan, Y. Sun, B. Zhou, C. Liu, D. Liu, W. Li and T. Wang, Nat. Energy, 2024, 9, 208–218.
15. R. Zhang, H. Chen, T. Wang, L. Kobera, L. He, Y. Huang, J. Ding, B. Zhang, A. Khasbaatar, S. Nanayakkara, J. Zheng, W. Chen, Y. Diao, S. Abbrent, J. Brus, A. H. Coffey, C. Zhu, H. Liu, X. Lu, Q. Jiang, V. Coropceanu, J.-L. Brédas, Y. Li, Y. Li and F. Gao, Nat. Energy, 2024, 10, 124–134.
16. L.-Y. Xu, W. Wang, X. Yang, S. Wang, Y. Shao, M. Chen, R. Sun and J. Min, Nat. Commun., 2024, 15, 1248.
17. Y. Jiang, S. Sun, R. Xu, F. Liu, X. Miao, G. Ran, K. Liu, Y. Yi, W. Zhang and X. Zhu, Nat. Energy, 2024, 9, 975–986.
18. Z. Zhao, S. Chung, Y. Y. Kim, M. Jeong, X. Li, J. Zhao, C. Zhu, S. Karuthedath, Y. Zhong, K. Cho and Z. Kan, Energy Environ. Sci., 2024, 17, 5666–5678.
19. N. Gao, P. Zhang, Z. Xu, S. Wu, L. Chen, J. Zhang, H. Bin and Y. Li, Small, 2025, 21, 2410679.
20. H. Liu, H. Bai, Y. Zhou, P. Li, W. Su, C. Liu, X. Liao, B. Song, X. Li, Z. Bi, C. Zhao, H. Liu, G. Lu, H. Du, L. Jiang, Y. Liu, R. Ma, W. Ma and Q. Fan, Mater. Sci. Eng., R, 2025, 162, 100879.
21. M. Liu, L. Wu, Y. Hai, Y. Luo, Y. Li, R. Chen, Y. Ma, T. Jia, Q. Li, S. Liu, R. Ma, Y.-P. Cai, J. Wu, G. Li and S. Liu, Adv. Mater., 2025, 37, 2503702.
22. J. Dong, Y. Li, C. Liao, X. Xu, L. Yu, R. Li and Q. Peng, Energy Environ. Sci., 2025, 18, 4982–4995.
23. M. Zhang, L. Zhu, J. Yan, X. Xue, Z. Wang, F. Eisner, G. Zhou, R. Zeng, L. Kan, L. Wu, W. Zhong, A. Zhang, F. Han, J. Song, N. Hartmann, Z. Zhou, H. Jing, H. Zhu, S. Xu, J. Nelson, Y. Zhang and F. Liu, Joule, 2025, 9, 101851.
24. Z. Zhao, S. Chung, L. Tan, J. Zhao, Y. Liu, X. Li, L. Bai, H. Lee, M. Jeong, K. Cho and Z. Kan, Energy Environ. Sci., 2025, 18, 2791–2803.
25. C. Li, J. Song, H. Lai, H. Zhang, R. Zhou, J. Xu, H. Huang, L. Liu, J. Gao, Y. Li, M. H. Jee, Z. Zheng, S. Liu, J. Yan, X.-K. Chen, Z. Tang, C. Zhang, H. Y. Woo, F. He, F. Gao, H. Yan and Y. Sun, Nat. Mater., 2025, 24, 433–443.
26. L. Liu, H. Li, J. Xie, Z. Yang, Y. Bai, M. Li, Z. Huang, K. Zhang and F. Huang, Adv. Mater., 2025, 37, 2500352.
27. H. Chen, Y. Huang, R. Zhang, H. Mou, J. Ding, J. Zhou, Z. Wang, H. Li, W. Chen, J. Zhu, Q. Cheng, H. Gu, X. Wu, T. Zhang, Y. Wang, H. Zhu, Z. Xie, F. Gao, Y. Li and Y. Li, Nat. Mater., 2025, 24, 444–453.
28. H. Fan, H. Yang, Y. Wu, C. Cui and Y. Li, Adv. Energy Mater., 2025, 15, 2405257.
29. L. Guo, J. Song, J. Deng, J. Qiao, J. Zhang, C. Li, S. Yuan, B. Han, M. H. Jee, Z. Ge, C. Zhang, G. Lu, X. Hao, H. Y. Woo and Y. Sun, Adv. Mater., 2025, 37, 2504396.
30. X. Dai, B. Fan, X. Xu and Q. Peng, Adv. Mater., 2025, 37, 2503072.
31. X. Xia, L. Mei, C. He, Z. Chen, N. Yao, M. Qin, R. Sun, Z. Zhang, Y. Pan, Y. Xiao, Y. Lin, J. Min, F. Zhang, H. Zhu, J.-L. Brédas, H. Chen, X.-K. Chen and X. Lu, J. Mater. Chem. A, 2023, 11, 21895–21907.
32. Z. Shen, J. Yu, G. Lu, K. Wu, Q. Wang, L. Bu, X. Liu, Y. Zhu and G. Lu, Energy Environ. Sci., 2023, 16, 2945–2956.
33. J. Wang, Y. Wang, M. Du, Y. Yu, C. Wang, W. Wang, Q. Guo, Y. Cui, S. Zhang and J. Hou, Energy Environ. Sci., 2024, 17, 8368–8378.
34. G. Zhang, X.-K. Chen, J. Xiao, P. C. Y. Chow, M. Ren, G. Kupgan, X. Jiao, C. C. S. Chan, X. Du, R. Xia, Z. Chen, J. Yuan, Y. Zhang, S. Zhang, Y. Liu, Y. Zou, H. Yan, K. S. Wong, V. Coropceanu, N. Li, C. J. Brabec, J.-L. Brédas, H.-L. Yip and Y. Cao, Nat. Commun., 2020, 11, 3943.
35. X. Liao, Q. Xie, Y. Guo, Q. He, Z. Chen, N. Yu, P. Zhu, Y. Cui, Z. Ma, X. Xu, H. Zhu and Y. Chen, Energy Environ. Sci., 2022, 15, 384–394.
36. X. Wu, X. Jiang, X. Li, J. Zhang, K. Ding, H. Zhuo, J. Guo, J. Li, L. Meng, H. Ade and Y. Li, Adv. Mater., 2023, 35, 2302946.
37. T. Wang, M. K. Ravva and J.-L. Brédas, Adv. Funct. Mater., 2016, 26, 5913–5921.
38. T. Wang and J.-L. Brédas, Matter, 2020, 2, 119–135.
39. T. Wang and J.-L. Brédas, J. Am. Chem. Soc., 2021, 143, 1822–1835.
40. R. A. Marcus, Rev. Mod. Phys., 1993, 65, 599.
41. V. Coropceanu, X.-K. Chen, T. Wang, Z. Zheng and J.-L. Brédas, Nat. Rev. Mater., 2019, 4, 689–707.
42. T. M. Burke, S. Sweetnam, K. Vandewal and M. D. McGehee, Adv. Energy Mater., 2015, 5, 1500123.
43. T. Wang, X. Chen, A. Ashokan, Z. Zheng, M. K. Ravva and J.-L. Brédas, Adv. Funct. Mater., 2018, 28, 1705868.
44. Z. Zheng, N. R. Tummala, T. Wang, V. Coropceanu and J.-L. Brédas, Adv. Energy Mater., 2019, 9, 1803926.
45. S. M. Pratik, G. Kupgan, J.-L. Brédas and V. Coropceanu, Energy Environ. Sci., 2025, 18, 841–852.
46. T. Unger, S. Wedler, F.-J. Kahle, U. Scherf, H. Bässler and A. Köhler, J. Phys. Chem. C, 2017, 121, 22739–22752.
47. M. Azzouzi, J. Yan, T. Kirchartz, K. Liu, J. Wang, H. Wu and J. Nelson, Phys. Rev. X, 2018, 8, 031055.
48. J. Jortner, J. Chem. Phys., 1976, 64, 4860–4867.
49. F.-J. Kahle, A. Rudnick, H. Bässler and A. Köhler, Mater. Horiz., 2018, 5, 837–848.
50. U. Rau, Phys. Rev. B: Condens. Matter Mater. Phys., 2007, 76, 085303.
51. B. Blank, T. Kirchartz, S. Lany and U. Rau, Phys. Rev. Appl., 2017, 8, 024032.
52. R. C. Hilborn, Am. J. Phys., 1982, 50, 982–986.
53. L. H. Smith, J. A. E. Wasey and W. L. Barnes, Appl. Phys. Lett., 2004, 84, 2986–2988.
54. R. Sun, J. Guo, Q. Wu, Z. Zhang, W. Yang, J. Guo, M. Shi, Y. Zhang, S. Kahmann, L. Ye, X. Jiao, M. A. Loi, Q. Shen, H. Ade, W. Tang, C. J. Brabec and J. Min, Energy Environ. Sci., 2019, 12, 3118–3132.
55. J. Wu, G. Li, J. Fang, X. Guo, L. Zhu, B. Guo, Y. Wang, G. Zhang, L. Arunagiri, F. Liu, H. Yan, M. Zhang and Y. Li, Nat. Commun., 2020, 11, 4612.
56. X. Li, A. Tang, H. Wang, Z. Wang, M. Du, Q. Guo, Q. Guo and E. Zhou, Angew. Chem., Int. Ed., 2023, 62, e202306847.
57. S. Plimpton, J. Comput. Phys., 1995, 117, 1–19.
58. W. L. Jorgensen and J. Tirado-Rives, J. Am. Chem. Soc., 1988, 110, 1657–1666.
59. N. E. Jackson, K. L. Kohlstedt, B. M. Savoie, M. Olvera De La Cruz, G. C. Schatz, L. X. Chen and M. A. Ratner, J. Am. Chem. Soc., 2015, 137, 6254–6262.
60. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, G. A. Petersson, H. Nakatsuji, X. Li, M. Caricato, A. V. Marenich, J. Bloino, B. G. Janesko, R. Gomperts, B. Mennucci, H. P. Hratchian, J. V. Ortiz, A. F. Izmaylov, J. L. Sonnenberg, D. Williams-Young, F. Ding, F. Lipparini, F. Egidi, J. Goings, B. Peng, A. Petrone, T. Henderson, D. Ranasinghe, V. G. Zakrzewski, J. Gao, N. Rega, G. Zheng, W. Liang, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, K. Throssell, J. A. Montgomery Jr, J. E. Peralta, F. Ogliaro, M. J. Bearpark, J. J. Heyd, E. N. Brothers, K. N. Kudin, V. N. Staroverov, T. A. Keith, R. Kobayashi, J. Normand, K. Raghavachari, A. P. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, J. M. Millam, M. Klene, C. Adamo, R. Cammi, J. W. Ochterski, R. L. Martin, K. Morokuma, O. Farkas, J. B. Foresman and D. J. Fox, Gaussian 16, Revision C.01, Gaussian, Inc., Wallingford CT, 2016.
61. M. J. Abraham, T. Murtola, R. Schulz, S. Páll, J. C. Smith, B. Hess and E. Lindahl, SoftwareX, 2015, 1–2, 19–25.
62. P. C. T. Souza, R. Alessandri, J. Barnoud, S. Thallmair, I. Faustino, F. Grünewald, I. Patmanidis, H. Abdizadeh, B. M. H. Bruininks, T. A. Wassenaar, P. C. Kroon, J. Melcr, V. Nieto, V. Corradi, H. M. Khan, J. Domański, M. Javanainen, H. Martinez-Seara, N. Reuter, R. B. Best, I. Vattulainen, L. Monticelli, X. Periole, D. P. Tieleman, A. H. De Vries and S. J. Marrink, Nat. Methods, 2021, 18, 382–388.
63. R. Alessandri, J. Barnoud, A. S. Gertsen, I. Patmanidis, A. H. De Vries, P. C. T. Souza and S. J. Marrink, Adv. Theory Simul., 2022, 5, 2100391.
64. J.-P. Ryckaert and A. Bellemans, Faraday Discuss. Chem. Soc., 1978, 66, 95.
65. R. Alessandri, J. J. Uusitalo, A. H. De Vries, R. W. A. Havenith and S. J. Marrink, J. Am. Chem. Soc., 2017, 139, 3697–3705.
66. S. F. Boys and F. Bernardi, Mol. Phys., 1970, 19, 553.
67. V. Coropceanu, J. Cornil, D. A. Da Silva Filho, Y. Olivier, R. Silbey and J.-L. Brédas, Chem. Rev., 2007, 107, 926–952.
68. B.-C. Lin, B. T. Koo, P. Clancy and C.-P. Hsu, J. Phys. Chem. C, 2014, 118, 23605–23613.
69. G. D'Avino, L. Muccioli, Y. Olivier and D. Beljonne, J. Phys. Chem. Lett., 2016, 7, 536–540.
70. B. Mennucci and J. Tomasi, J. Chem. Phys., 1997, 106, 5151–5158.
71. R. J. Cave and M. D. Newton, Chem. Phys. Lett., 1996, 249, 15–19.
72. T. Wang, V. Coropceanu and J.-L. Brédas, Chem. Mater., 2019, 31, 6239–6248.

---

Footnote

† These authors contribute equally to this work.

---