Organic solar cells: evolution of morphological and electronic properties as a function of steric hindrance from the acceptor side groups

Abstract

High non-radiative recombination losses are currently limiting the efficiency of organic solar cells (OSCs). Attaching sterically bulky side groups to acceptors is a strategy that has been recently developed to reduce such losses. However, it remains unclear how these side groups impact the nanoscale morphology and electronic properties in neat acceptor films and donor:acceptor blends as a function of their steric hindrance due to the experimental challenges in accessing these aspects, which hinders the further development of more efficient acceptors. Here, we characterize these aspects using a tight combination of all-atom molecular dynamics simulations and long-range corrected density functional theory calculations. As representative systems, we select PB2 as the polymer donor and three Y6-core-based acceptors functionalized with benzene, trimethylbenzene, and triisopropylbenzene groups (denoted as BTP-B, BTP-Bme, and BTP-Biso, respectively). As steric hindrance increases, excessive acceptor aggregation is suppressed; acceptor_acceptor packing distances increase, reducing electronic couplings and intermolecular electron-transfer rates (k_e). Interestingly, the PB2_PB2 packing distances also increase; however, the packing pattern changes, enhancing electronic couplings and the interchain hole-transfer rate (k_h). This results in a more balanced k_e/k_h ratio in the PB2:BTP-Bme blend. Additionally, the PB2_acceptor packing distances increase, which elevates the charge-transfer states and decreases electronic couplings between these states and the ground state, thereby lowering the corresponding non-radiative recombination rates and losses. Together, these yield the lowest non-radiative recombination losses in the PB2:BTP-Bme-based OSCs. Overall, we present a comprehensive picture that describes how such a strategy improves the nanoscale morphology and electronic properties of OSCs, which can be beneficial for designing highly efficient acceptors.

---

Broader context

Compared with inorganic and perovskite solar cells, organic solar cells (OSCs) suffer from a higher (∼0.15 eV) non-radiative energy loss, which largely limits their open-circuit voltage (V_OC) and ultimately their efficiency. Attaching sterically bulky side groups to acceptors is a strategy that has been recently developed to reduce these losses via increasing [decreasing] the radiative [non-radiative] recombination rates. Since accessing the nanoscale morphology is experimentally challenging, it remains poorly understood how these sterically bulky side groups impact the nanoscale morphology and the related electronic properties of neat acceptor films and donor:acceptor blends. For solving this issue a robust computational methodology that combines long-range corrected density functional theory calculations and all-atom molecular dynamics simulations is used. Using this multiscale method, we draw a comprehensive picture that describes how such a strategy improves the nanoscale morphology and electronic properties of OSCs, which can be beneficial for designing highly efficient acceptors.

Introduction

Organic solar cells (OSCs) have become a prominent photovoltaic technology, which combines advantages such as light weight, semi-transparency, mechanical flexibility, and large-scale processability.^1–3 Although the emergence and development of non-fullerene small-molecule acceptors (NF-SMAs) have pushed the OSC power conversion efficiencies (PCEs) to values approaching 21%,^4–15 these PCEs still lag behind those of inorganic and perovskite solar cells largely due to the high non-radiative recombination losses (ΔE_nr) typical of OSCs.^5,6,13–19

ΔE_nr can be quantitatively determined by the following expression:^20,21

where k_B denotes the Boltzmann constant, T is the device temperature, EQE_EL is the electroluminescence quantum efficiency, p_e indicates the probability of a photon, generated through a radiative recombination process, escaping from the device, and k_r and k_nr are the radiative and non-radiative recombination rates, respectively. This relationship clearly points out the expectation that minimizing ΔE_nr requires maximizing k_r and minimizing k_nr. To increase the k_r values, attaching sterically bulky side groups to acceptors, such as triphenylamine,¹⁴ tetra-phenylethylene,¹⁵ benzene,^18,22,23 triisopropyl-benzene,^18,22 2Br-CzP,²⁴ acridine,²⁵ and norbornadiene groups,^26–28 has recently been demonstrated to be an effective strategy.

Regarding the evaluation of k_nr from the charge-transfer (CT) states to the ground state, the Marcus–Levich–Jortner model is often used:^29–31

where ħ represents the reduced Planck's constant, V_el is the electronic coupling between a CT state and the ground state, S_qm denotes the Huang–Rhys parameter with S_qm = λ_qm/ħω_qm, λ_qm represents the quantum-mechanical contribution to the reorganization energy described by means of an effective vibrational mode ħω_qm, λ_c is the classical contribution to the reorganization energy, and E^ad_CT denotes the adiabatic (relaxed) CT energy. Accordingly, a smaller k_nr calls for a lower reorganization energy, a lower electronic coupling, and a higher CT-state energy.

To reduce ΔE_nr by increasing k_r and decreasing k_nr, in a recent work, Hou and co-workers incorporated benzene, trimethylbenzene, and triisopropylbenzene side groups to the shoulder positions of the Y6 core moiety and synthesized the BTP-B, BTP-Bme, and BTP-Biso acceptors, respectively (Fig. 1).¹⁸ These side groups were able to improve the photoluminescence quantum yield (PLQY) and thus k_r, as their bulkiness hindered the torsion between the end and core moieties of the Y6 acceptor; this improves the acceptor planarity and rigidity and thus effectively reduces the reorganization energy. In going from BTP-B to BTP-Bme and then to BTP-Biso, the ΔE_nr values of the PB2:acceptor-based OSCs evolved from 0.253 eV to 0.190 eV and to 0.229 eV, and the corresponding PCEs from 6.26% to 18.5% and to 15.4%.^18,32 The steric hindrances brought by the trimethylbenzene groups into BTP-Bme have been found experimentally to assist the formation of a favorable global morphology (i.e., optimal degree of phase separation and domain size), which gives the PB2:BTP-Bme-based OSCs the highest short-circuit current density (J_SC), fill factor (FF), and thus PCE.¹⁸ In addition, introducing BTP-Bme as a third component into the PB2:PBDB-TF binary blend also yielded a high efficiency of 19.3%.¹⁸

Fig. 1 Chemical structures of PB2, BTP-B, BTP-Bme, and BTP-Biso.

Although a previous work¹⁸ has provided reference systems to understand how bulky groups on the sides of NF-SMAs can alter the active-layer morphology, electronic properties, and device performance as a function of steric hindrance, many aspects remain unexplored, owing to experimental limitations in accessing morphological and electronic properties at the nanoscale. For instance, in the neat acceptor films: (i) how do the intramolecular conformation and the corresponding site energy (i.e., electron affinity, EA) evolve; (ii) how do the intermolecular packing (i.e., packing pattern, distance, and order); and the related electron-transfer rates vary; and in the PB2:acceptor blend films: (iii) how do the intermolecular packing (i.e., acceptor_acceptor, PB2_PB2, and PB2_acceptor packing) and the relevant electronic properties change (i.e., the electron-transfer rates between acceptor molecules, hole-transfer rates between PB2 chains, and the interfacial CT states and their radiative/non-radiative recombination processes) as a function of the degree of steric hindrance brought into the acceptors. Here, to address these aspects, we combine long-range corrected density functional theory (DFT) calculations and molecular dynamics (MD) simulations to characterize the morphological and electronic properties of pure acceptor systems (BTP-B, BTP-Bme, and BTP-Biso) and the PB2:acceptor blend systems (PB2:BTP-B, PB2:BTP-Bme, and PB2:BTP-Biso) at the nanoscale. Doing so allows us to draw a comprehensive picture that describes how such a strategy (i.e., introducing side groups with appropriate steric hindrance into NF-SMAs) improves the morphological and electronic properties of the OSC active layers, which can help the further design of efficient NF-SMA materials.

Results and discussion

Intramolecular conformations, intermolecular packing configurations, and related electronic properties in neat acceptor films

As a first step, we examined the impact of steric hindrance induced by side groups on the intramolecular conformations in the neat BTP-B, BTP-Bme, and BTP-Biso films. Generally, the intramolecular conformation of an acceptor molecule that consists of the end and core moieties (Fig. S1a) depends on: (i) the rotations between the end moieties and the core moiety since both are composed of rigid and coplanar fused rings, and (ii) the intermolecular interactions between a given molecule and the surrounding molecules. To this end, we first evaluated the evolution of molecular potential energies by scanning the dihedral angle that characterizes the rotation between the acceptor end and core moieties (i.e., S–C–C–C as highlighted in Fig. S1b).³³ As shown in Fig. 2a, for all the BTP-B, BTP-Bme, and BTP-Biso molecules, the cis conformation possesses the lowest potential energy. As steric hindrance increases from BTP-B to BTP-Bme and then to BTP-Biso, it requires more energy to rotate the end moiety with respect to the core moiety around the global minimum. In addition, the energy barrier to the transition from the cis to trans conformation becomes larger. From the perspective of intramolecular interactions, it can be anticipated that in going from BTP-B to BTP-Bme and then to BTP-Biso, the acceptor molecule becomes more rigid and coplanar, and has a greater probability of adopting the cis conformation.

Fig. 2 (a) Evolution of potential energies by rotating the end moiety with respect to the core moiety; (b) representative end_end, end_core, and core_core π–π packing configurations extracted from the MD-simulated BTP-B film; (c) radial distribution functions for the end_end, end_core, and core_core interactions and (d) end_end, end_core and core_core π–π packing fractions; and (e) distributions of the S–C–C–C dihedral angles between the end and core moieties in the MD-simulated BTP-B, BTP-Bme, and BTP-Biso films.

We now focus on intermolecular interactions. Fig. S2 presents the radial distribution functions (RDFs) for the acceptor backbone–backbone interactions in the MD-simulated BTP-B, BTP-Bme, and BTP-Biso films; here, the RDF measures how the relative density of atoms varies as a function of distance away from a reference atom [if ρ_global is the average number density of atoms, the local averaged density at a distance r is ρ_local = ρ_global·g(r)].^34–36 As the steric hindrance increases from BTP-B to BTP-Bme and then to BTP-Biso, the first g(r) peak that relates to the nearest-neighbor intermolecular packing weakens and shifts to a larger distance. The weakening of the first g(r) peak suggests a reduction in the π–π packing density. To confirm, we extracted and obtained 280, 216, and 125 pairs with π–π packing features from the MD-simulated BTP-B, BTP-Bme and BTP-Biso films, which corresponds to a decrease in π–π packing density from 0.451 nm⁻³ to 0.318 nm⁻³ and then to 0.158 nm⁻³. This evolution is expected to inhibit any excessive aggregation of acceptor molecules and thus to lead to a smaller degree of phase separation in the donor:acceptor blends. The shift of the first g(r) peak towards a larger distance then indicates an increase in the π–π packing distance; this result is in line with the Grazing Incidence Wide-Angle X-ray Scattering (GIWAXS) measurements that show that the (010) d-spacing in the neat BTP-B, BTP-Bme, and BTP-Biso films rises from 3.59 Å to 3.65 Å and then to 3.87 Å.¹⁸ These variations are consistent with the fact that increasing steric hindrance decreases the intermolecular interaction energies. Indeed, as summarized in Table S1, the intermolecular backbone–backbone interaction energies in the MD-simulated films decreased from −17.2 ± 6.2 to −13.8 ± 4.4 and then to −13.0 ± 3.8 kcal mol⁻¹ as the steric hindrance increases from BTP-B to BTP-Bme and then to BTP-Biso.

Considering that the acceptor backbone is composed of the end and core moieties, there may exist three types of intermolecular π–π packing configurations in neat films, i.e., end_end (end moieties on top of one another), end_core (end moiety on top of the core moiety), and core_core (core moieties on top of one another)^37–40 (Fig. 2b). Accordingly, Fig. 2c presents the partial RDFs for the end_end, end_core, and core_core interactions in the MD-simulated BTP-B, BTP-Bme, and BTP-Biso films. Since the steric hindrance induced by the bulky side groups on the core moieties prevents any significant intermolecular π–π packing around the core moieties, it is expected that the preferential packing follows the same order in the three films: end_end > end_core > core_core, which is indeed reflected by the relative g(r) heights around ∼4 Å. In addition, the first g(r) peaks relevant for the end_core and core_core interactions diminish markedly with increased steric hindrance. To quantitatively characterize these findings, Fig. 2d provides the fractions of each packing configuration in the three films. In all the cases, the fractions of these packing configurations follow the sequence end_end > end_core > core_core. In addition, as steric hindrance increases, the proportions of the end_core and core_core packing configurations markedly decrease; especially, the latter decreases from 9% to 1%, then to 0. Accordingly, the end_end fraction increases substantially and even surpasses 90% in the BTP-Biso film. As displayed in Fig. 2b, in the end_end packing configuration, the two end moieties that interact with each other do not cross over the S–C–C–C dihedral angle between the end and core moieties, and thus impact only very weakly the rigidity and planarity of the whole molecule.

Based on the above discussion, it is anticipated that: (i) both the cis and trans conformations could be present in neat BTP-B, BTP-Bme, and BTP-Biso films, with cis being the dominant conformation; and (ii) as the steric hindrance goes up, the proportion of cis molecules will augment. To verify these points, we analyzed the distributions of the S–C–C–C dihedral angles that characterize the torsion between the acceptor core and end moieties,³⁵ in the MD-simulated BTP-B, BTP-Bme, and BTP-Biso films. As shown in Fig. 2e, the BTP-B, BTP-Bme, and BTP-Biso molecules can adopt both the cis and trans conformations (cis and trans correspond to the distributions at around 0–30° and 130–170°, respectively, with the former being more planar than the latter). As expected, with an increase in the degree of steric hindrance, the proportion of cis conformations greatly augments, which results in narrower and more uniform distributions for intramolecular conformations and more rigid and coplanar acceptor backbones.

As a result, the energetic disorder of the transport levels (i.e., EAs, for the acceptor molecules) is expected to decrease in going from BTP-B to BTP-Bme and then to BTP-Biso; we recall that a small energetic disorder often plays a positive role in determining the electronic processes occurring in OSCs.^23,41–43 To better assess this aspect, we extracted 900 BTP-B, 900 BTP-Bme, and 900 BTP-Biso molecules from the MD simulations and computed their EA values. Gaussian fits to the probability distributions of the EA values provide the total energetic disorder σ_T (Fig. 3a). Indeed, as the steric hindrance goes up from BTP-B to BTP-Bme and then BTP-Biso, the σ_T value decreases from 57.4 meV to 51.8 meV and then to 49.9 meV (Fig. 3b); we note that all these values are smaller than that found in PC₇₁BM (i.e., σ_T = 77.0 meV).^41,44

Fig. 3 (a) Total disorders for EAs computed from the distribution in EA values of 900 molecules; the average EA value, <EA>, was used as a reference for each system. (b) Total disorders (σ_T), dynamic disorders (σ_D), and static disorders (σ_S) of EAs in the MD-simulated BTP-B, BTP-Bme, and BTP-Biso films.

To explore more deeply the origin of the evolution of σ_T values as a function of steric hindrance, we divided σ_T into its dynamic (σ_D) and static (σ_S) components via the relationship of σ_T² = σ_D² + σ_S² (see more details in the Computational Methodology).^41,45,46 The former is related to the electron–vibration interactions that lead to a time-dependent variation of EAs with a standard deviation σ_D; the latter is associated with the variations in the positions of the molecules (or the lack of perfect order of the molecules) in their neat films, which results in a time-independent modulation of EAs with a standard deviation σ_S. In line with the increase in the extent of planarity and rigidity from BTP-B to BTP-Bme and BTP-Biso, the values of σ_D and σ_S both decrease; however, the latter decreases more significantly, which indicates that the smaller σ_T values in BTP-Bme and BTP-Biso stem mainly from a remarkable decrease in their static disorder.

Having in hand the 280, 216, and 125 π–π packing pairs extracted from the MD-simulated BTP-B, BTP-Bme, and BTP-Biso films, we evaluated the electron transfer rates (k_e) between adjacent acceptor molecules. In the framework of semi-classical Marcus theory, the k_e values are jointly determined by the electronic coupling (or transfer integral, V_e), reorganization energy (λ_e), and energetic mismatches between EAs (ΔEA). A larger k_e value calls for a higher V_e, a lower λ_e, as well as a |ΔEA| value approaching λ_e. As collected in Table S3, as the steric hindrance increases from BTP-B to BTP-Bme and BTP-Biso, there occurs a decrease in the λ_e value from 294 meV to 289 meV and 286 meV, a change in |ΔEA| from 68 meV to 64 meV and to 66 meV, and a decrease in the V_e value from 12.7 meV to 10.1 meV and 9.7 meV. Overall, their combined effects result in a decrease in the k_e value from 2.38 × 10¹² s⁻¹ to 1.42 × 10¹² s⁻¹ and then to 1.19 × 10¹² s⁻¹ (for transfer from an acceptor molecule with a smaller EA to one with a larger EA) and from 2.02 × 10¹¹ s⁻¹ to 1.57 × 10¹¹ s⁻¹ and then to 1.19 × 10¹¹ s⁻¹ (for transfer from an acceptor with a larger EA to one with a smaller EA). Table S4 summarizes the fractions of the π–π packing pairs corresponding to various orders of magnitude in these rates. At most of the higher levels of electron transfer rates (i.e., ≥10¹³, ≥10¹², ≥10¹¹, and ≥10¹⁰), the fractions of π–π packing pairs consistently follow the sequence of BTP-B > BTP-Bme > BTP-Biso. These results are consistent with the experimental findings that the electron mobility (µ_e) decreases from 5.27 × 10⁻³ cm² V⁻¹ s⁻¹ for BTP-B to 3.31 × 10⁻³ cm² V⁻¹ s⁻¹ for BTP-Bme and to 1.91 × 10⁻³ cm² V⁻¹ s⁻¹ for BTP-Biso.¹⁸

Acceptor_acceptor packing configurations and related electronic properties in the PB2:acceptor blend films

At this stage, we wish to examine whether the above findings for neat acceptor films can be translated to blend films in the presence of the PB2 polymer donor. Fig. 4a collects the RDFs for the acceptor_acceptor backbone interactions in the MD-simulated PB2:BTP-B, PB2:BTP-Bme, and PB2:BTP-Biso blends. Here as well, the acceptor π–π packing distance [density] is found to increase [decrease] with steric hindrance, as reflected by the right shift and weakening of the first g(r) peak. From the MD-simulated PB2:BTP-B, PB2:BTP-Bme, and PB2:BTP-Biso blends, we can extract 143, 101, and 53 acceptor π–π packing pairs, respectively, which further confirms a decrease in the π–π packing density from 0.138 nm⁻³ to 0.090 nm⁻³ and to 0.042 nm⁻³.

Fig. 4 (a and b) Radial distribution functions for the acceptor_acceptor and PB2_PB2 backbone interactions in the MD-simulated PB2:BTP-B, PB2:BTP-Bme and PB2:BTP-Biso blends; (c) illustrations of the PB2 D and A moieties; and (d) partial radial distribution functions for the PB2 D_D, D_A, and A_A interactions in the MD-simulated PB2:BTP-B, PB2:BTP-Bme, and PB2:BTP-Biso blends.

Fig. S3a–c display the partial RDFs for the acceptor end_end, end_core, and core_core interactions in the MD-simulated PB2:BTP-B, PB2:BTP-Bme, and PB2:BTP-Biso blends, respectively. The preferential packing follows exactly the same order as that in the neat acceptor films, i.e., end_end > end_core > core_core, as reflected by comparing the heights of the relevant g(r) peaks around ∼4 Å. In addition, the first g(r) peaks related to the interactions involving the core moieties are again markedly suppressed as the steric hindrance increases. Fig. S3d provides the quantitative fractions for each packing configuration in the three blends, which follow the sequence end_end > end_core > core_core; a larger steric hindrance results in the proportions of end_core and core_core packings decreasing significantly and a marked increase in the proportion of the end_end packing configuration (coming close to 0.90). Overall, the effects of steric hindrance on intermolecular packing configurations in the neat acceptor films appear to be largely transferable to the blend films.

Table S5 lists the intermolecular electron-transfer rates, k_e, between two adjacent acceptors in the blends as well as the related parameters. With steric hindrance going up from BTP-B to BTP-Bme and BTP-Biso, the averaged k_e values decrease from 1.43 × 10¹² s⁻¹ to 1.30 × 10¹² s⁻¹ and 1.12 × 10¹² s⁻¹ (for transfer from an acceptor molecule with a smaller EA to one with a larger EA), and from 1.62 × 10¹¹ s⁻¹ to 1.39 × 10¹¹ s⁻¹ and 1.22 × 10¹¹ s⁻¹ (for transfer from an acceptor molecule with a larger EA to one with a smaller EA); this evolution is mainly attributed to a decrease in electronic couplings. In addition, directional electron transport pathways are present in the MD-simulated PB2:BTP-B and PB2:BTP-Bme blends (Fig. S4), which further supports their efficient local electron transfer properties. The µ_e values were measured experimentally to be 1.03 × 10⁻⁵ cm² V⁻¹ s⁻¹, 1.15 × 10⁻⁴ cm² V⁻¹ s⁻¹, and 6.21 × 10⁻⁵ cm² V⁻¹ s⁻¹ for the PB2:BTP-B, PB2:BTP-Bme, and PB2:BTP-Biso blend films, respectively. We recall that the µ_e value mainly depends on two factors: local intermolecular electron transfer rates and global morphology.⁴⁸ Therefore, the discrepancy in the sequences of the evaluated electron transfer rates (i.e., BTP-B > BTP-Bme > BTP-Biso) and the measured µ_e values (i.e., BTP-Bme > BTP-Biso > BTP-B) is expected to be related to the unfavorable global morphology of the PB2:BTP-B film. Indeed, atomic force microscopy (AFM) measurements have revealed that the PB2:BTP-B film surface is covered by large granular clusters of acceptor aggregates, which leads to a much larger value of the root mean square surface roughness (R_q) of the PB2:BTP-B film vs. those of the PB2:BTP-Bme and PB2:BTP-Biso films (i.e., 90.8 nm vs. 1.01 nm and 0.92 nm).¹⁸ This excessive aggregation of BTP-B molecules results in inefficient exciton dissociation and significant charge recombination and thus negatively impacts electron mobility.¹⁸

PB2_PB2 packing configurations and related electronic properties in the PB2:acceptor blend films

Here as well, as the steric hindrance increases, the PB2_PB2 interchain packing distance becomes larger, as reflected by the right shift of the first g(r) peak from 4.45 Å to 4.55 Å and then to 4.65 Å in the RDFs for the PB2_PB2 backbone interactions in the MD-simulated PB2:BTP-B, PB2:BTP-Bme and PB2:BTP-Biso blends (see Fig. 4b). However, a major difference with respect to what we discovered regarding the acceptor_acceptor packing is that, with the increase in the acceptor steric hindrance, the PB2_PB2 packing intensity shows a pronounced enhancement. The reason is that the increase in acceptor steric hindrance can not only weaken the PB2_acceptor interactions but also create additional free volume, which provides a greater possibility of PB2_PB2 interactions. Indeed, we were able to extract 48, 52, and 57 PB2_PB2 pairs with π–π packing features from the MD-simulated PB2:BTP-B, PB2:BTP-Bme, and PB2:BTP-Biso blends, respectively.

Based on the PB2_PB2 π–π packing pairs extracted from the MD-simulated blends, Table S6 displays the interchain hole transfer rates, k_h, and their related parameters in the framework of semi-classical Marcus theory.⁴⁷ Going from BTP-B to BTP-Bme and to BTP-Biso slightly increases the k_h value from 1.02 × 10¹² s⁻¹ to 1.31 × 10¹² s⁻¹ and to 1.87 × 10¹² s⁻¹ (for transfer from a PB2 chain with a smaller ionization potential, IP, to one with a larger IP), and from 0.80 × 10¹¹ s⁻¹ to 1.15 × 10¹¹ s⁻¹ and to 1.62 × 10¹¹ s⁻¹ (for transfer from a PB2 segment with a larger IP to one with a smaller IP). Table S7 summarizes the fractions of the PB2_PB2 π–π packing pairs corresponding to different orders of magnitude in k_h values. At most of the higher levels of hole transfer rates (i.e., ≥10¹³, ≥10¹², ≥10¹¹, and ≥10¹⁰), the fractions of π–π packing pairs consistently follow the sequence PB2:BTP-Biso > PB2:BTP-Bme > PB2:BTP-B. This correlates with the increase in electronic coupling between adjacent PB2 chains, from 12.4 meV to 15.4 meV and to 15.8 meV (Table S6). In contrast to the finding that the increase in acceptor_acceptor π–π packing distance decreases the electronic coupling for the electron transfer process, the increase in the PB2_PB2 π–π packing distance instead increases the electronic coupling for the hole transfer process. Actually, the electronic coupling depends not only on the intermolecular packing distance but also on the intermolecular relative displacements along the backbones.^49,50 In fact, in going from PB2:BTP-B to PB2:BTP-Bme and then to PB2:BTP-Biso, the PB2_PB2 packing pattern (i.e., the moieties, D or A, through which the PB2 chains pack, Fig. 4c) changes, namely from D_D (D moieties on top of one another) > D_A (D moiety on top of A moiety) > A_A (A moieties on top of one another) to D_D ≈ A_A > D_A and then to D_D > A_A > D_A, as shown in Fig. 4d. While the calculated interchain hole transfer rate follows the sequence PB2:BTP-Biso > PB2:BTP-Bme > PB2:BTP-B, the hole mobility (µ_h) measured experimentally is as follows: PB2:BTP-Bme > PB2:BTP-Biso > PB2:BTP-B (i.e., 1.28 × 10⁻⁴ cm² V⁻¹ s⁻¹ > 1.04 × 10⁻⁴ cm² V⁻¹ s⁻¹ > 9.95 × 10⁻⁵ cm² V⁻¹ s⁻¹). This difference can be attributed to the lack of sufficient phase-separated morphology of the PB2:BTP-Biso active layer. This feature is supported by the much smaller Flory–Huggins interaction parameter, χ, computed for the PB2:BTP-Biso blend (0.026 K) than for the PB2:BTP-B (0.120 K) and PB2:BTP-Bme blends (0.053 K), which is detrimental to exciton dissociation and charge transport⁵¹ (we note that these χ values are calculated via the expression , where K is a positive constant and γ_D/γ_A denotes the surface energy of the neat donor/acceptor film).^18,39,52

Among the three blends, PB2:BTP-Bme presents the k_e/k_h ratio closest to 1 (i.e., 0.99 vs. 1.40 and 0.60 for transfer from smaller to larger EA/IP and 1.21 vs. 2.03 and 0.75 for transfer from larger to smaller EA/IP for the PB2:BTP-Bme blend vs. the PB2:BTP-B and PB2:BTP-Biso blends). The combination of this appropriate k_e/k_h ratio with the favorable active-layer morphology leads to a more balanced µ_e/µ_h of 0.90 in the PB2:BTP-Bme-based OSCs,¹⁸ which decreases the charge recombination and is thus conducive to increasing both the J_SC (i.e., 25.1 mA cm⁻² for PB2:BTP-Bme vs. 12.9 mA cm⁻² for PB2:BTP-B and 21.4 mA cm⁻² for PB2:BTP-Biso) and FF (i.e., 0.799 for PB2:BTP-Bme vs. 0.580 for PB2:BTP-B and 0.762 for PB2:BTP-Biso).¹⁸

PB2_acceptor packing configurations and related electronic properties in the PB2:acceptor blend films

The increase in the acceptor steric hindrance is also found to enhance the PB2_acceptor packing distance, as seen from the right shift of the first g(r) peak position from 4.35 Å to 4.55 Å and 4.75 Å (Fig. 5a) in the RDFs of the PB2_acceptor backbone interactions in the MD-simulated PB2:BTP-B, PB2:BTP-Bme and PB2:BTP-Biso blends. This larger intermolecular distance is expected to weaken the electrostatic interaction between the excited electron and hole, thereby increasing the CT-state energy.^46,53 To address this point, we evaluated the energies of the lowest singlet CT electronic states for the PB2_acceptor π–π packing pairs extracted from the MD-simulated PB2:BTP-B, PB2:BTP-Bme, and PB2:BTP-Biso blends. Fig. 5b displays the normalized energetic distributions of the CT states for the three blends, which are fitted with the Gaussian functions. Indeed, in going from PB2:BTP-B to PB2:BTP-Bme and PB2:BTP-Biso, the average energies of the CT states, E_CT^avg, increase from 1.455 to 1.504 and 1.521 eV, and their standard deviations, σ, decrease from 0.122 to 0.109 and 0.107 eV.

Fig. 5 (a) Radial distribution functions for the PB2_acceptor backbone interactions; (b) normalized energetic distributions of the CT states fitted with Gaussian functions; (c) non-radiative recombination rates k_nr; and (d) non-radiative recombination losses ΔE_nr from the CT states to the ground state for the MD-simulated PB2:BTP-B, PB2:BTP-Bme, and PB2:BTP-Biso blends.

Based on the Marcus–Levich–Jortner model that takes quantum effects into consideration,^29–31 Fig. 5c displays the non-radiative recombination rates, k_nr, from the lowest singlet CT states to the ground state, where the quantum-mechanical component of the reorganization energy, λ_qm, was set to be close to half of the intramolecular reorganization energy (i.e., 0.11 eV). The k_nr values follow the order PB2:BTP-B > PB2:BTP-Bme > PB2:BTP-Biso; this can be assigned both to the increase in the CT-state energies and to the decrease in electronic couplings (related to the enhanced PB2_acceptor packing distance) and reorganization energy values (correlating with the increased acceptor planarity and rigidity), with the former playing a more prominent role,^54,55 i.e., PB2:BTP-B < PB2:BTP-Bme < PB2:BTP-Biso, see the parameters listed in Table S8.

Furthermore, we evaluated the radiative recombination rates, k_r, from the CT states to the ground state through the Einstein coefficient relation.²¹ As shown in Table S9, the PB2:BTP-Bme and PB2:BTP-Biso blends have slightly higher k_r values than the PB2:BTP-B blend (i.e., 4.19 × 10⁵ s⁻¹ and 3.78 × 10⁵ s⁻¹ vs. 2.60 × 10⁵ s⁻¹). This comes from the higher CT-state energies and transition dipole moments between the CT states and the ground state for the PB2:BTP-Bme and PB2:BTP-Biso blends.

Based on the relationship given in the Introduction and the calculated k_nr and k_r values, the non-radiative recombination losses, ΔE_nr, relevant for processes from the CT states to the ground state were evaluated and illustrated in Fig. 5d for the MD-simulated PB2:BTP-B, PB2:BTP-Bme, and PB2:BTP-Biso blends (we note that the p_e value is taken to be 0.25, a value that lies in the range of the largest values of ∼0.2–0.3 measured for the light-outcoupling coefficient in organic light-emitting devices).^35,56 The ΔE_nr values follow the order of PB2:BTP-B > PB2:BTP-Bme > PB2:BTP-Biso, which is in line with the sequence in experimental V_OC values: PB2:BTP-B (0.837 V) < PB2:BTP-Bme (0.927 V) < PB2:BTP-Biso (0.945 V).¹⁸ We note, however, that the experimental ΔE_nr value for the PB2:BTP-Bme (0.190 eV) blend is smaller than those of the PB2:BTP-B (0.253 eV) and PB2:BTP-Biso blends (0.229 eV).¹⁸ This can be attributed to other beneficial aspects of the PB2:BTP-Bme blend, such as higher µ_e and µ_h values, a more balanced µ_e/µ_h ratio, and a more favorable active-layer morphology, which further reduce the ΔE_nr value in the PB2:BTP-Bme-based OSC.

Finally, we correlated the donor_acceptor packing pattern with the k_nr and ΔE_nr values in the MD-simulated PB2:BTP-Bme blend. Considering that the PB2 backbone consists of D and A moieties (Fig. 4c) and BTP-Bme is composed of end and core moieties (Fig. S1a), there exist four possible donor_acceptor packing configurations, i.e., D_end, D_core, A_end, and A_core, as illustrated in Fig. 6a. Fig. 6b underlines that the different types of PB2_BTP-Bme π–π packing configurations in fact possess different k_nr values, i.e., A_core > D_end > D_core > A_end. Accordingly, the A_end packing configuration has the lowest ΔE_nr value while the A_core packing configuration carries the highest ΔE_nr value (Fig. 6c). Therefore, efforts to increase [decrease] the proportion of A_end [A_core] packing configurations would be beneficial to further reduce non-radiative recombination losses in the PB2:BTP-Bme system. In addition, we evaluated the exciton-dissociation rates for the D_end, D_core, A_end, and A_core packing patterns (Table S10), and found that there also exist differences among the four packing patterns (i.e., 2.88 × 10⁹ s⁻¹ for D_end, 1.34 × 10⁹ s⁻¹ for D_core, 2.07 × 10⁹ s⁻¹ for A_end, and 1.37 × 10¹¹ s⁻¹ for A_core, respectively).

Fig. 6 (a) Representative D_end, D_core, A_end, and A_core π–π packing configurations extracted from the MD-simulated PB2:BTP-Bme blend; (b) non-radiative recombination rates k_nr; and (c) non-radiative recombination losses ΔE_nr from the CT states to the ground state for the four types of packing configurations in the MD-simulated PB2:BTP-Bme blend.

Conclusions

Examining the morphological and electronic properties in neat acceptor films and PB2:acceptor blends, via a combination of long-range-corrected DFT calculations and MD simulations, has enabled us to draw a comprehensive picture that describes the molecular-scale effects related to acceptor side groups as a function of the steric hindrance they induce. The main results are as follows:

(i) In neat acceptor films, as the steric hindrance increases from BTP-B to BTP-Bme and BTP-Biso, the intramolecular torsion potentials become steeper and the intermolecular interactions weakened, which leads to more rigid and coplanar acceptor backbones and thus decreased EA energetic disorders. The overall intermolecular packing order is effectively enhanced by suppressing acceptor core-related packing patterns. The decrease in intermolecular interactions not only impairs acceptor aggregation, which avoids oversized acceptor aggregates, but also enlarges the π–π packing distance without altering the packing pattern, which decreases electronic couplings and thus intermolecular electron transfer rates.

(ii) In PB2:acceptor blends, the steric-hindrance effects on the acceptor are similar to those found in neat acceptor films. Although the PB2 interchain packing distance increases along with the increase in acceptor steric hindrance, the PB2 interchain packing pattern does change, which accounts for an increase in electronic couplings and thus interchain hole transfer rates. Overall, these effects result in the most balanced k_e/k_h ratio in the PB2:BTP-Bme blend, which contributes to more balanced µ_e/µ_h and higher J_SC and FF.

(iii) The PB2_acceptor packing distance also becomes larger, which increases the CT-state energies and lowers the electronic couplings related to the non-radiative recombination processes from the CT states to the ground state. In conjunction with the decreased reorganization energy, these aspects reduce the non-radiative recombination rates and the related energy losses from the CT states to the ground state. Combining the favorable characteristics of morphological and electronic properties leads to the PB2:BTP-Bme blend displaying the lowest overall ΔE_nr value; the A_end packing configuration in this blend is also found to provide the lowest non-radiative recombination rate and related energy loss.

Overall, the methodology that combines long-range corrected DFT calculations and all-atom MD simulations can also be used to examine the local morphology and related electronic properties in other polymer/NF-SMA blends. In addition, these results are expected to translate well to other polymer/NF-SMA systems in which the sterically bulky side groups are attached to NF-SMAs to reduce non-radiative recombination losses. However, due to the temporal and spatial limitations of all-atom MD simulations, the methodology cannot be directly used to quantitatively describe the global morphology and its long-term evolution. To this end, the coarse-grained MD is an appropriate choice, as it reduces the degree of freedom of the system and extends the temporal and spatial scales close to experimental conditions.

Computational methodology

Molecular dynamics (MD) simulations

All the MD simulations were carried out with the LAMMPS package⁵⁷ and the force field of optimized potentials for liquid simulations-all atom (OPLS-AA).^34,35,58,59 To accurately capture the intramolecular and intermolecular interactions for the PB2, BTP-B, BTP-Bme, and BTP-Biso components, the force field was reparametrized using long-range corrected DFT calculations: (i) the atomic partial charges were derived by fitting the electrostatic potential (ESP) evaluated at the ωB97XD/cc-PVTZ level of theory; (ii) the bond lengths and angles were directly taken from the geometries optimized at the ωB97XD/6-31G(d,p) level of theory while the harmonic force constants were kept as they are; and (iii) the parameters for the inter-ring dihedrals along the PB2 chains as well as the BTP-B, BTP-Bme, and BTP-Biso molecules were fitted in light of their torsion potentials computed at the ωB97XD/6-31G(d,p) level of theory. Here, the related 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.

To validate the reparametrized force field, we carried out a direct comparison between the DFT- and force field-calculated (from the fitted parameters) torsion potentials and a good agreement was found (including the minimum, maximum, and profile), as shown in Fig. S5.

The initial models for neat films were built by randomly placing 300 BTP-B, 300 BTP-Bme, and 300 BTP-Biso molecules in three cubic cells with a low density of 0.02 g cm⁻³, respectively. Similarly, the initial models for the blend films were built by randomly placing 41 PB2 chains and 242 BTP-B molecules, 44 PB2 chains and 242 BTP-Bme molecules, and 48 PB2 chains and 242 BTP-Biso molecules in three cubic boxes with a low density of 0.02 g cm⁻³, with each PB2 chain consisting of 8 repeat units. The choices of the PB2 chain length, number of PB2 chains, and number of BTP-B/BTP-Bme/BTP-Biso molecules were based on three considerations: (i) the experimental PB2:acceptor weight ratio of 1 : 1;¹⁸ (ii) the computational feasibility as well as the temporal and spatial limitations of MD simulations, and (iii) the chain length of 5–8 repeat units were also successfully used to examine the local morphology and electronic properties in other polymer:NF-SMA blends.^39,61 The MD simulations were first performed with the NPT (constant number of molecules, pressure, and temperature) ensemble for 30 ns at a temperature of 600 K and a pressure of 1 atm. Then, all neat and blend films were rapidly cooled from 600 K to 400 K, with a cooling rate of ∼50 K ns⁻¹. To simulate the experimental annealing process, the following MD simulations were carried out: (i) at 373.15 K for 20 ns, where 373.15 K is the experimental annealing temperature;¹⁸ (ii) at 350 K for 10 ns; and (iii) at 298.15 K for 20 ns. A cutoff of 12 Å was used for the summation of van der Waals interactions and the particle–particle particle-mesh (PPPM) solver was used for 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.

Density functional theory (DFT) calculations

The intermolecular backbone–backbone interaction energies for the acceptor_acceptor π–π packing pairs were calculated at the ωB97XD/6-31G(d,p) level of theory, with the geometries for all pairs kept as those extracted from the MD-simulated neat acceptor films and all side chains substituted with methyl groups. To avoid overstabilization of these energies induced by basis set superposition error (BSSE), the counterpoise correction method developed by Boys and Bernardi was adopted.⁶²

In order to examine the energetic disorders for the EA of BTP-B, BTP-Bme, and BTP-Biso in their neat films, three frames were extracted at 0.5 ns intervals from the last 1.5 ns for each MD-simulated acceptor film. To obtain the total disorder (σ_T), the EA values for 900 BTP-B/BTP-Bme/BTP-Biso molecular structures (300 molecules × 3 frames) were first evaluated; the σ_T value was then obtained by a Gaussian fitting to the EA-energy distribution. The dynamic disorder (σ_D) was obtained by considering the geometries of BTP-B, BTP-Bme, and BTP-Biso molecules generated from the MD trajectories as a function of time. The coordinates of 200 structures for each of 8 randomly selected molecules were extracted every 14 fs from the MD simulation trajectories. Following this, the σ_D values were obtained through a Gaussian fit to the EA-energy distribution for each of the 8 randomly selected molecules. The values summarized in Fig. 3b correspond to the average of the 8 randomly selected molecules. The static disorders (σ_S) are then simply calculated as:.^41,46

Ideally, to accurately capture the electrostatic and polarization effects from neighboring molecules on the EA values, the quantum mechanics/molecular mechanics (QM/MM) approach is preferred. However, given that the acceptor molecules contain at least 165 atoms and the MD-simulated films include more than 49 500 atoms, such an approach would be computationally prohibitive. Therefore, in this context, we evaluated the EA values by combining single-molecule DFT calculations with the polarizable continuum model (PCM), which may underestimate the EA value but is not expected to change the trend among the BTP-B, BTP-Bme, and BTP-Biso molecules.

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

where ħ denotes the reduced Planck's constant; V_i, the electronic coupling between the initial state and final state; λ_i, the reorganization energy; and ΔE_i, the energy difference between the initial state and final state. V_i between the initial and final states was taken as the coupling between the frontier molecular orbitals (i.e., the lowest unoccupied molecular orbital for the electron and the highest occupied molecular orbital for the hole) of the two molecules in the π–π packing pair; V_i was evaluated with the fragment orbital method.⁵⁰ ΔE_i between the initial and final states was estimated as the energy difference between the EAs or ionization potential (IP) (i.e., EA for electron and IP for hole) of the two molecules in the π–π packing pair. The λ_i reorganization energy is composed of intra- (λ_intra) and inter-molecular (λ_inter) components. 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),^63,64 considering that there are currently no straightforward models to accurately estimate the λ_inter value in solid-state environments. This value of 0.1 eV was also widely adopted for studying analogous OSC systems, which has been found to reasonably reproduce the experimental trends.^40,65 All the DFT calculations were carried out 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 while considering a typical dielectric constant of 3.5).^66,67

The Marcus–Levich–Jortner model was adopted to evaluate the non-radiative recombination rates, k_nr from the CT states to the ground state:^29–31

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 calculated in the framework of the generalized Mulliken–Hush approach.⁶⁸ Here, the E^ad_CT value was evaluated as: E^ad_CT = E^avg_CT − λ, where E^avg_CT is the averaged vertical transition energies of the CT states computed for the PB2_acceptor π–π packing pairs at their ground-state geometries as extracted from the 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.^63,64 Considering that the PB2_acceptor π–π packing pairs extracted from the 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⁻¹.^35,69 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 adopted the λ_qm value of 0.11 eV (a typical value close to half the value of λ_intra); then, λ_c was taken as λ_c = λ − λ_qm. Here, all the 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 ε₀ denotes the vacuum permittivity, c is the speed of light and µ_CT→S0, the transition dipole moment between the CT state and the ground state.

Author contributions

Ailing Yin: writing – original draft, validation, methodology, formal analysis. Ziwen Yu: validation. Jingyi Ren: methodology. Jiacong Feng:formal analysis. 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

Additional details relevant to the molecular structures, RDFs and electronic properties discussed in the main text are provided in the supplementary information (SI). Supplementary information is available. See DOI: https://doi.org/10.1039/d6el00013d.

Acknowledgements

The work at Jilin University was financially supported by the National Key R&D Program of China (No. 2023YFB3003000), 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.
2. G. Yu, J. Gao, J. C. Hummelen, F. Wudl and A. J. Heeger, Science, 1995, 270, 1789.
3. 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.
4. 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.
5. 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.
6. J. Dong, Y. Li, C. Liao, X. Xu, L. Yu, R. Li and Q. Peng, Energy Environ. Sci., 2025, 18, 4982.
7. K. Hu, Y. Ge, H. Yang, Y. Xu, J. Qian, X. Zhu, Y. Wu, C. Cui and Y. Li, Energy Environ. Sci., 2025, 18, 9194.
8. 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.
9. H. Mou, Y. Yin, H. Chen, J. Xu, J. Ding, C. Ju, J. Zhu, Y. Wang, W. Chen, G. Xu, T. Zhang, J. Li, Y. Li and Y. Li, J. Am. Chem. Soc., 2025, 147, 21241.
10. L. Wang, C. Chen, Z. Gan, J. Cheng, Y. Sun, J. Zhou, W. Xia, D. Liu, W. Li and T. Wang, Adv. Mater., 2025, 37, 2419923.
11. Y. Xie, J. Tian, X. Yang, J. Chen, S. Yu, D. Tang, X. Hu, Y. Sun and M. Lv, Adv. Mater., 2025, 37, e02485.
12. X. Duan, J. Zhang, J. Kong, B. Song, J. Qiao, M. H. Jee, J. Deng, W. Li, H. Y. Woo, X. Hao, G. Lu, J. Song and Y. Sun, Adv. Mater., 2025, 38, e14076.
13. L. Guo, L. Wu, T. Jia, H. Zhang, J. Song, X. Xie, M. H. Jee, H. Ma, S. Liu, G. Lu, H. Y. Woo, Z. Wang, F. Gao and Y. Sun, Angew. Chem., Int. Ed., 2025, 64, e202516421.
14. Y. Chen, X. Duan, J. Zhang, Z. Ge, H. Ma, X. Sun, H. Zhang, J. Gao, X. Wang, X. Wang, Z. Tang, R. Yang, F. Gao and Y. Sun, Energy Environ. Sci., 2025, 18, 6214.
15. J. Wang, T. Liu, Y. Cheng, C. Sun, X. Jiang, J. Zhao, Y. Li, S. Ying, Y. Liu, H. Jiang, H. Lu, Y. Liu, Z. Bo, X. Bao and S. Yan, Adv. Funct. Mater., 2026, 36, e24391.
16. I. Ramirez, M. Causa, Y. Zhong, N. Banerji and M. Riede, Adv. Energy Mater., 2018, 8, 1703551.
17. Q. Liu and K. Vandewal, Adv. Mater., 2023, 35, 2302452.
18. J. Ren, S. Zhang, Z. Chen, T. Zhang, J. Qiao, J. Wang, L. Ma, Y. Xiao, Z. Li, J. Wang, X. Hao and J. Hou, Angew. Chem., Int. Ed., 2024, 63, e202406153.
19. Y. Lin, Z. Ma and Z. Tang, Mater. Horiz., 2025, 12, 7714.
20. U. Rau, Phys. Rev. B, 2007, 76, 085303.
21. B. Blank, T. Kirchartz, S. Lany and U. Rau, Phys. Rev. Appl., 2017, 8, 024032.
22. L. Ma, S. Zhang, J. Zhu, J. Wang, J. Ren, J. Zhang and J. Hou, Nat. Commun., 2021, 12, 5093.
23. Y. Shi, L. Zhu, Y. Yan, M. Xie, G. Liang, J. Qiao, J. Zhang, X. Hao, K. Lu and Z. Wei, Adv. Energy Mater., 2023, 13, 2300458.
24. Y. Pan, L. Guo, M. H. Jee, G. Dai, Z. Ge, J. Zhang, X. Duan, J. Song, X. Li, H. Y. Woo and Y. Sun, Adv. Energy Mater., 2024, 16, 2403747.
25. J. Liu, X. Duan, J. Zhang, Z. Ge, L. Liu, J. Qiao, Y. Li, Z. Bi, H. Zhang, J. Gao, J. Yan, S. Liu, J. Zhang, Z. Tang, X. Hao, W. Ma, F. Gao and Y. Sun, Angew. Chem., Int. Ed., 2025, 64, e202500129.
26. H. Lu, W. Liu, H. Jin, H. Huang, Z. Tang and Z. Bo, Adv. Funct. Mater., 2022, 32, 2107756.
27. C. e. Zhang, R. Zheng, H. Huang, G. Ran, W. Liu, Q. Chen, B. Wu, H. Wang, Z. Luo, W. Zhang, W. Ma, Z. Bo and C. Yang, Adv. Energy Mater., 2024, 14, 2303756.
28. X. Zhu, C. Gu, Y. Cheng, H. Lu, X. Wang, G. Ran, W. Zhang, Z. Tang, Z. Bo and Y. Liu, Adv. Mater., 2025, 37, 2507529.
29. J. Jortner, J. Chem. Phys., 1976, 64, 4860.
30. T. Unger, S. Wedler, F.-J. Kahle, U. Scherf, H. Bässler and A. Köhler, J. Phys. Chem. C, 2017, 121, 22739.
31. F.-J. Kahle, A. Rudnick, H. Bässler and A. Köhler, Mater. Horiz., 2018, 5, 837.
32. T. Zhang, C. An, Y. Cui, J. Zhang, P. Bi, C. Yang, S. Zhang and J. Hou, Adv. Mater., 2022, 34, 2105803.
33. T. Wang and J.-L. Brédas, Adv. Funct. Mater., 2019, 29, 1806845.
34. T. Wang, M. K. Ravva and J.-L. Brédas, Adv. Funct. Mater., 2016, 26, 5913.
35. T. Wang and J.-L. Brédas, Matter, 2020, 2, 119.
36. T. Wang and J.-L. Brédas, J. Am. Chem. Soc., 2021, 143, 1822.
37. G. Han, T. Hu and Y. Yi, Adv. Mater., 2020, 32, 2000975.
38. G. Kupgan, X.-K. Chen and J.-L. Brédas, Mater. Today Adv., 2021, 11, 100154.
39. Y. Huang, Z. Yu, A. Yin, G. F. Han, X. Y. Lang, J.-L. Brédas, T. Wang and Q. Jiang, Angew. Chem., Int. Ed., 2024, 64, e202418225.
40. Z. Yu, T. Zhang, Y. Li, G.-F. Han, X.-Y. Lang, J.-L. Brédas, T. Wang and Q. Jiang, EES Sol., 2025, 1, 1017.
41. G. Kupgan, X.-K. Chen and J.-L. Brédas, ACS Mater. Lett., 2019, 1, 350.
42. N. A. Ran, J. A. Love, C. J. Takacs, A. Sadhanala, J. K. Beavers, S. D. Collins, Y. Huang, M. Wang, R. H. Friend, G. C. Bazan and T.-Q. Nguyen, Adv. Mater., 2016, 28, 1482.
43. Y. L. Lin, M. A. Fusella and B. P. Rand, Adv. Energy Mater., 2018, 8, 1702816.
44. S. Xie, Y. Xia, Z. Zheng, X. Zhang, J. Yuan, H. Zhou and Y. Zhang, Adv. Funct. Mater., 2018, 28, 1705659.
45. N. R. Tummala, Z. Zheng, S. G. Aziz, V. Coropceanu and J.-L. Brédas, J. Phys. Chem. Lett., 2015, 6, 3657.
46. Z. Zheng, N. R. Tummala, T. Wang, V. Coropceanu and J.-L. Brédas, Adv. Energy Mater., 2019, 9, 1803926.
47. R. A. Marcus, Rev. Mod. Phys., 1993, 65, 599.
48. T. M. Clarke and J. R. Durrant, Chem. Rev., 2010, 110, 6736.
49. J.-L. Brédas, J. P. Calbert, D. A. da Silva Filho and J. Cornil, Proc. Natl. Acad. Sci. U.S.A., 2002, 99, 5804.
50. V. Coropceanu, J. Cornil, D. A. da Silva Filho, Y. Olivier, R. Silbey and J.-L. Brédas, Chem. Rev., 2007, 107, 926.
51. Y. Cai, Q. Li, G. Lu, H. S. Ryu, Y. Li, H. Jin, Z. Chen, Z. Tang, G. Lu, X. Hao, H. Y. Woo, C. Zhang and Y. Sun, Nat. Commun., 2022, 13, 2369.
52. 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.
53. J. Wang, X. Jiang, H. Wu, G. Feng, H. Wu, J. Li, Y. Yi, X. Feng, Z. Ma, W. Li, K. Vandewal and Z. Tang, Nat. Commun., 2021, 12, 6679.
54. X.-K. Chen, V. Coropceanu and J.-L. Brédas, Nat. Commun., 2018, 9, 5295.
55. V. Coropceanu, X.-K. Chen, T. Wang, Z. Zheng and J.-L. Brédas, Nat. Rev. Mater., 2019, 4, 689.
56. L. H. Smith, J. A. E. Wasey and W. L. Barnes, Appl. Phys. Lett., 2004, 84, 2986.
57. S. Plimpton, J. Comput. Phys., 1995, 117, 1.
58. W. L. Jorgensen and J. Tirado-Rives, J. Am. Chem. Soc., 1988, 110, 1657.
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.
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. J. Zhou, C. Guo, L. Wang, C. Chen, Z. Gan, Y. Sun, C. Liu, J. Zhou, Z. Chen, D. Gao, W. Xia, D. Liu, T. Wang and W. Li, Nano Energy, 2024, 129, 109988.
62. S. F. Boys and F. Bernardi, Mol. Phys., 1970, 19, 553.
63. B.-C. Lin, B. T. Koo, P. Clancy and C.-P. Hsu, J. Phys. Chem. C, 2014, 118, 23605.
64. G. D'Avino, L. Muccioli, Y. Olivier and D. Beljonne, J. Phys. Chem. Lett., 2016, 7, 536.
65. S. M. Pratik, G. Kupgan, J.-L. Brédas and V. Coropceanu, Energy Environ. Sci., 2025, 18, 841.
66. B. Mennucci and J. Tomasi, J. Chem. Phys., 1997, 106, 5151.
67. T. Wang, X.-K. Chen, A. Ashokan, Z. Zheng, M. K. Ravva and J.-L. Brédas, Adv. Funct. Mater., 2018, 28, 1705868.
68. R. J. Cave and M. D. Newton, Chem. Phys. Lett., 1996, 249, 15.
69. T. Wang, V. Coropceanu and J.-L. Brédas, Chem. Mater., 2019, 31, 6239.

---