Towards first-principles prediction of open-circuit voltage and short-circuit current density in small-molecule BHJ solar cells

Abstract

We present a theoretical framework to predict the open-circuit voltage and short-circuit current density of small-molecule bulk heterojunction solar cells. The method is based on the reciprocity relation and on a quantum mechanical description of the elementary processes. Radiative and non-radiative decay rates are evaluated through Fermi's golden rule. The Franck–Condon weighted density of states, computed via the generating function formalism, is employed as a genuine predictive tool, in that it enters the expression for both radiative and non-radiative decay rates and determines the absorption profile governing photogeneration. The method requires only a limited set of experimentally accessible inputs, namely the redox potentials and absorption spectra of the donor and acceptor materials, and involves a small number of adjustable parameters with negligible impact on the results. This makes the framework suitable for studying different donor–acceptor combinations in small-molecule organic solar cells based on either fullerene or non-fullerene acceptors. The model provides a direct interpretation of open-circuit voltage by linking voltage losses and photocurrent generation to specific processes, such as non-radiative recombination, subgap effects, and absorption-edge broadening. It also captures the temperature dependence of open-circuit voltage and short-circuit current density, which is essential for evaluating device performance under realistic operating conditions. Overall, the proposed approach provides a unified and computationally efficient framework for predicting open-circuit voltage and short-circuit current density and for relating them to the underlying physical processes in small-molecule organic solar cells.

---

Introduction

Organic solar cells (OSCs) offer compelling advantages over traditional inorganic photovoltaic technologies,^1,2 including tunable optoelectronic properties,^3,4 mechanical flexibility,^5,6 lightweight construction, compatibility with low-temperature solution-processing,⁷ and the potential for scalable, cost-effective manufacturing. Such attributes make OSCs promising candidates for emerging applications such as wearable electronics,^8,9 building-integrated photovoltaics, and semi-transparent power-generating surfaces.¹⁰

Among the most explored OSC architectures, the bulk heterojunction (BHJ) design represents a significant improvement over earlier bilayer structures, which suffer from limited exciton diffusion and inefficient charge dissociation.^1,11 In BHJ devices, donor and acceptor materials are intimately mixed on the nanoscale, allowing for improved exciton harvesting and charge separation efficiency.^12–14

Although polymer donor-small molecule acceptor systems represent an important class of high-efficiency organic solar cells,¹⁵ bulk heterojunction devices entirely based on small molecules have achieved comparable power conversion efficiencies (PCE), exceeding 14% and, in some cases, approaching or surpassing 17%.^16–20 Unlike polymers, small molecules enable more consistent control over the active layer morphology in BHJs. This improved morphological control enhances charge transport and reduces recombination losses, and also contributes to better device reproducibility.^21,22

Despite their promising characteristics, organic solar cells have not yet achieved commercial success, due to challenges related to stability and scalability, but primarily because of performance limitations such as lower power conversion efficiencies compared to inorganic solar cells,^23–25 even though a steady increase in the maximum PCE has been observed over the years, breaking lately the 20% threshold.²⁶

A solar cell is characterized in terms of open-circuit voltage, V_OC, and short-circuit current density, J_SC. The factors that determine these key quantities have been the subject of extensive investigation and are now well established.^27–34 Unlike their inorganic counterparts, organic materials suffer from low dielectric constants, leading to strongly bound excitons whose dissociation requires energetic offsets at the donor–acceptor interface, which inherently limits V_OC.³⁵ Additional losses arise from trap-assisted recombination,³⁶ non-radiative decay,^14,30 and subgap absorption associated with charge-transfer (CT) states and tail states, which introduce recombination pathways below the optical gap.^37,38 Moreover, the carbon-based backbone of organic materials supports high-frequency vibrational modes, primarily C–C stretching vibrations,³⁹ which can couple with electronic transitions and promote non-radiative decay, especially from CT states. These vibronic effects, together with limited charge carrier mobility and interfacial energy mismatches, reduce the overall efficiency of the device.^28,40

To overcome this issue, considerable effort has been devoted to identify new materials for solar cells.^2,30 Machine learning algorithms are becoming increasingly popular.^41–45 Despite their effectiveness, these methods rely heavily on extensive experimental data, which can still limit their applicability.^46,47 Therefore, developing traditional approaches alongside machine learning is important, as it offers additional insight into device operation and supports improved OSC design.^27,28,48

With this goal in mind, we present a theoretical approach aimed at (1) predicting the open-circuit voltage and short-circuit current of bulk heterojunction solar cells based on non-polymeric donor and acceptor materials, and (2) interpreting the main factors that govern the open-circuit voltage. The method computes, from first principles, the most relevant quantities needed to V_OC and J_SC.

Our approach is novel in its use of the Franck–Condon weighted density of states (FCWD), which plays a central role in this work. While this quantity has previously been used to rationalize voltage loss mechanisms within simplified frameworks,^27,28 we extend its use to a genuinely predictive level. Specifically, we predict the emission and non-radiative decay rates of the charge-transfer state, together with the entire absorption profile, by computing the FCWD from first principles, without introducing any empirical parameter.

Furthermore, the Franck–Condon weighted density of states inherently incorporates thermal effects,^49,50 which enables a direct investigation of the temperature dependence of V_OC and J_SC. This capability is particularly relevant for assessing solar-cell performance under realistic outdoor operating conditions.^51,52

All the quantities needed to predict V_OC and J_SC are herein obtained by quantum chemical computations with the following exceptions: (i) the charge-transfer state energy is evaluated from experimental redox free energies estimated by voltammetry and/or photoelectron spectroscopy; (ii) the donor and acceptor excitation energies are inferred from the observed absorption spectra of the isolated components, and assigned by means of the computed FCWD. Notably, measuring redox potentials and absorption spectra on isolated components is far simpler than fabricating the full solar cell. Finally, the method uses four adjustable parameters: the number density of active CT complexes, layer thickness, refractive index, and blend density, but the results are rather insensitive to their values.

In order to test the method, we applied it to eight devices for which the experimentally reported blend composition, donor:acceptor ratio, V_OC, and J_SC are available.

While most calculations presented here have been performed at a fixed temperature of 300 K, since experimental values of the open-circuit voltage and short-circuit current density are available only at room temperature for the investigated systems, we have also carried out calculations at different temperatures to quantitatively assess the ability of our method to reproduce the temperature dependence of V_OC and J_SC.

Theoretical methodology

When the solar cell is at open-circuit, the overall current flow through the device is null. Hence, recombination and photocurrent densities must coincide:⁵³

J_rec(V_OC) = J_ph(V_OC), (1)

where J_rec(V) is due to the migration of electrons and holes that, ultimately, recombine at the interface, whereas J_ph(V) is the difference between the dark and light current at the internal voltage V. In order to determine V_OC, we adopt the Shockley non-ideal diode equation:^54,55where T is the absolute temperature, q the elementary charge, and k_B the Boltzmann constant. J₀ is the saturation current density, i.e., the sum of the current densities arising from all radiative and non-radiative recombination mechanisms operating within the device:

J₀ = J_0,rad + J_0,nr. (3)

The subscript 0 is used because, even in the dark, photons are absorbed as the solar cell remains in equilibrium with the thermal radiation from the surroundings.^27,48 Finally, n_id is the ideality factor; that empirical quantity equals 1 when only direct radiative recombination occurs, it is expected to range between 1 and 2 for trap-assisted recombination, and has been reported to be lower than 1 when recombination takes place at the interface between the absorber layer and the electrode in thin-film cells.^56,57

At open circuit, the internal voltage V is equal to V_OC. Substituting eqn (1) into (2) and solving for V_OC, we obtain:

As is well known, the actual V_OC can be decomposed into a radiative component and a term accounting for non-radiative losses:^27,58

V_OC = V_OC,rad − ΔV_OC,nr. (5)

For an ideal device where only radiative processes are present, we have J₀ = J_0,rad and n_id = 1. According to eqn (4):²⁷where J_ph(V_OC,rad) is the photocurrent when only radiative processes are considered. The quantity V_OC,rad represents the maximum possible value of V_OC in a system where radiative recombination is the only active mechanism. Note that V_OC,rad is lower than the Shockley–Queisser limit unless the external quantum efficiency (EQE), i.e., the ratio of the number of charge carriers collected at the electrodes by the solar cell to the number of incident photons at a given frequency is a step function of the radiation energy E:
with Q^(e) being external quantum efficiency and E_g the optical bandgap.⁵³

When the superposition principle applies,^59,60 i.e., the total current under illumination is the sum of photocurrent and dark current, it follows that

J_ph(V_OC) = J_SC, (7)

where J_SC is the short circuit current density, which becomes

J_ph(V_OC,rad) = J_SC (8)

under the conditions of eqn (6). Therefore, we can use the equality (8) in eqn (6):We further note that V_OC (and a fortiori V_OC,rad) must be much higher than k_BT/q for the device to operate with reasonable photovoltaic efficiency. By applying this condition to eqn (6), and taking into account eqn (8), V_OC,rad can be safely approximated asJ_SC is typically measured with greater accuracy than J_ph(V_OC), as it corresponds to the current under conditions of negligible resistance. Indeed, under open-circuit, the absence of external current prevents the series resistance from affecting the device behavior. Even when the conditions for the superposition principle are not strictly fulfilled, approximation in eqn (7) introduces only a minor error in the calculation of V_OC, of the order of k_BT/q.⁵³

The short circuit current density, in the Shockley–Queisser approach, is defined as

where ϕ_sun(E) is the solar emission spectrum.⁶¹

The present model is simplified in that it assumes infinite carrier mobility and neglects bulk transport phenomena such as bimolecular recombination, trap-assisted Shockley–Read–Hall losses, and field-dependent charge extraction. Consequently, the calculated J_SC is expected to be somewhat overestimated, at least for D:A blends with limited carrier mobility or significant transport losses.

Under thermal equilibrium, emission originates from carriers at temperature T, resulting in a spectral shape close to blackbody radiation. Therefore, the radiative saturation current density is:^53,61

where ϕ_BB(E) is the spectral photon flux obtained by dividing Planck's blackbody radiation by photon energy:

Assuming perfect electron transport and no losses of photoinduced charge carriers, the absorptance of the device, A(E) is equivalent to its ideal EQE when operating as a solar cell. Therefore we can write:

Q^(e)(E) = A(E) (14)

in eqn (11) and (12). In light of eqn (9)–(14), it is essential to obtain accurate expressions for ϕ_sun and A(E) in order to compute V_OC,rad and J_SC. We adopt the AM1.5G standard reference spectrum as ϕ_sun.^62,63 Assuming no reflectance at the surface, unity reflectance from the back electrode and invoking Lambert–Beer law with no interference, we can express the absorptance as:

A(E) = 1 − e^−2α(E)d, (15)

where d is the thickness of the layer, herein taken as an adjustable parameter with the default value of 100 nm, α the absorption coefficient and the factor 2 accounts for the double pass of the beam due to complete reflection. Although we make use of eqn (15) to compute A(E), absorptance can be approximated as

A(E) ≈ 2α(E)d (16)

for very thin films.

Following previous work, we express the absorption coefficient α(E) as:^27,64

α(E) = α_A(E) + α_D(E) + α_CT(E), (17)

where the suffixes A, D, and CT refer to the absorption of donor, acceptor and CT complex at interface, respectively.

In most of the spectrum, the absorption of donor:acceptor blends can be approximated as a linear combination of the individual spectra of the donor and acceptor, weighted by their blending ratios. However, this approximation can fail in the lower-energy region, where the absorption of both D and A can be vanishingly small, yet absorption shoulders or low intensity peaks are possibly observed. These localized tail states within the bandgap are commonly attributed to structural defects and imperfections inherent to the low crystallinity characteristic of thin-film bulk heterojunction solar cells. These defects can lead to local fluctuations in the donor's HOMO–LUMO energy levels, potentially giving rise to trap states. Alternatively, such features may originate from absorption processes associated with charged polarons, phenomena invoked in cells utilizing polymeric donors with extended conjugation in combination with fullerene acceptors.^38,64 However, one must also consider the contribution of tails associated with intermolecular charge-transfer absorption, which involves direct excitation from the donor's HOMO to the acceptor's LUMO, as illustrated in Fig. 1a. These tails are detected when there is a large energy offset between the absorption edges of the donor and/or the acceptor and the CT state energy. Otherwise, the CT absorption is buried under the much stronger signals of the individual components.⁶⁵

Fig. 1 (a) Scheme of the sub band-gap process, i.e., the intermolecular charge-transfer transition from the HOMO of the donor to the LUMO of the acceptor. (b) Portion of the Jablonski diagram presenting the formation of a charge-transfer state via the photo-electron transfer (PET) process from the excited donor molecule, which could be followed by either radiative or non radiative charge recombination to the ground state. The corresponding radiative and non-radiative decay rates of the CT state are denoted by k_rad and k_nr, respectively.

Although the absorption cross sections of both CT states and impurity-induced transitions are noticeably smaller than those of singlet excitons, their contribution to photovoltaic losses is significant and cannot be omitted because it can be detrimental for power conversion efficiency.^27,38

The absorption coefficient of the species S (D or A) in a blend can be described under the Born–Oppenheimer and Condon approximations. The derivation also invokes the relationship between absorption and emission as described by the Einstein coefficients, assumes the validity of the Lambert–Beer law, adopts the electric dipole approximation,⁶⁶ and includes the isotropic orientational averaging of transition dipole moments. The resulting expression for α_S is given in SI units as:

Here, ω is the angular frequency of the radiation, w_S is the weight fraction of the species S, M_S is its molar mass, is Avogadro's number, n is the refractive index (real part, assumed constant), ρ is the blend density, μ^S_tr is the transition electric dipole moment for absorption leading to an excited state with electronic energy E_S above the ground state, and Ω_S(E_S − ħω, T) is the FCWD for the absorption.

In our model, the refractive index n and the density ρ are treated as adjustable parameters. For the former, we use the default value n = 1.5.²⁷ The blend density is set per default to 1.2 kg dm⁻³, because ρ ∼ 1.0 kg dm⁻³ for BHJ cells, with reported values ranging from 0.8 to 1.5 kg dm⁻³.⁶⁷

For donor and acceptor species, the weight fraction w_S is derived from the actual composition of a given device. The dipole strength (|μ^S_tr|)² is estimated at the time-dependent density functional theory (TDDFT) level (see the computational details) and corresponds to the absorption leading to the first singlet excited state for the systems considered here. However, our methodology is also capable of treating cases where multiple excited states, including triplet states, are involved.^14,68–70

The Franck–Condon weighted density of vibrational states is computed using the generating function approach (GF), which accounts for the entire vibrational manifold, includes thermal effects, and incorporates mode mixing (see the computational details).^31,69–74 Equilibrium geometries and Hessian matrices of the involved electronic states are obtained by DFT/TDDFT computations. Even E_S could in principle be obtained by TDDFT calculations. However, for the sake of accuracy, we extract excitation energies from absorption spectra of pure donor or acceptor films, by aligning the computed and observed Franck–Condon weighted densities, as illustrated below.^31,75

We apply eqn (18) also to the absorption leading to the charge-transfer state, with the following details: (1) the molar mass of the active CT complex is taken as the sum of the molar masses of the donor and the acceptor. (2) Its weight fraction is set such that the number density of the active CT complex is 100 times lower than the minimum number density of the donor and acceptor, consistent with previous studies assuming for the CT exciton a volume 100 times larger than that of a localized exciton.²⁷ The weight fractions of D and A are then adjusted so that the sum of all three weight fractions equals one. (3) The transition dipole moment, μ^CT_tr, connecting the ground and CT states (Fig. 1a) is evaluated at the TDDFT level by using the diabatization scheme detailed in the computational details. Calculations are performed for the minimum-energy configuration of the D:A complex in the ground state, as determined by molecular mechanics computations. (4) Under the assumption of weak interaction between D and A, the Franck–Condon weighted density of states of the CT state corresponds to that of a bimolecular process, and is thus obtained by convolving the FCWDs of the two half-reactions in eqn (19):

The convolution is given by

where Ω_D/D⁺ and Ω_A/A⁻ are the FCWDs for the processes in eqn (19), whose global reaction energy equals the CT energy (Fig. 1):

E_CT = IP(D) − EA(A), (21)

in which IP and EA denote ionization potential and electron affinity, respectively. Both quantities are here inferred from available experimental data. In summary, eqn (6)–(21) allow the computation of V_OC,rad and J_SC.

Non-radiative voltage losses, represented by the term ΔV_OC,nr (eqn (5)), are better recognized when formulated in terms of reciprocity relations that link the device operating as a solar cell (under illuminated conditions) to its electroluminescent behavior as a light-emitting diode (LED) in the dark. The light-emitting diode external quantum efficiency Q^(e)_LED of the solar cell is defined as the ratio between the number of photons emitted outside the device by the equivalent LED (when operating in the dark) and the number of CT states generated within it by the application of a voltage bias. This quantity can be expressed in terms of the recombination current densities:

where J_rad(V) is the current density due to charge-recombination events that produce photons, while J_nr(V) represents the current density due to non-radiative recombination events that do not emit photons but lead to dissipation at the bias voltage V.

Considering a bias voltage V = V_OC and using the Shockley law for J_rad together with the approximation used in eqn (10):

Inserting eqn (23) in eqn (22):Using (24), (1) and (7):The latter equation, in conjunction with (10) and (5), gives:Eqn (26) and (10) illustrate the effects of the reciprocity between absorption and emission: A high-performance solar cell combines efficient photon absorption to enhance J_SC and hence V_OC,rad with effective radiative emission, as the latter mitigates the effects of nonradiative losses and thus enhances V_OC.⁴⁰

Q^(e)_LED(V_OC) can be related to the internal quantum efficiency for luminescence, Q⁽ⁱ⁾_LED(V_OC), which is defined as the ratio between the number of photons generated inside the LED and the number of CT states formed in the device under the application of the voltage bias V_OC. By disregarding surface recombination and assuming that recombination rates are independent of the spatial location within the absorber – an assumption consistent with the use of an isotropic and spatially uniform refractive index – Q⁽ⁱ⁾_LED can be expressed in terms of radiative (k_rad) and non-radiative (k_nr) recombination rates of the CT state, as schematically illustrated in Fig. 1b:^61,76

Eqn (27) and (22) can be related. Both J_nr and k_nr quantify photon losses through non-radiative processes, i.e., they account for recombination events that do not result in light emission. In contrast, only a fraction, denoted p_e, of the photons emitted via radiative recombination – and thus contributing to k_rad – are externally emitted, thereby also contributing to J_rad. This consideration allows us to write:

This expression can be readily understood by interpreting p_e as the probability that a photon generated by radiative recombination escapes the device without being reabsorbed. Inserting (28) into (26) yields:²⁷Upon neglecting interference effects and the spatial dependence of rates, p_e turns out to be:^61,77In the thin-film limit, it follows from eqn (16) and (30) that

Radiative recombination rates are herein computed using Fermi's golden rule, with the FCWD obtained by the GF approach:^31,71

Ω^em_CT is the FCWD for the emission of the CT complex giving rise to charge recombination, computed through the convolution of the FCWDs for the inverse half-reactions in eqn (19). Here, the transition dipole moment for the vertical emission of CT, (μ^CT_tr), should be evaluated for the equilibrium geometry that the active D:A complex assumes in the excited state. Nevertheless, the same transition dipole moment as for absorption is herein used, since determining the excited-state equilibrium geometry of the D:A complex at the TDDFT level would be too demanding for such large systems. However, this approximation should not introduce a major error, as in the solid phase the D:A pair is expected to undergo only minimal structural relaxation upon excitation.⁷⁸

Within Fermi's golden rule formalism, the non-radiative recombination rate turns out to be:

where V_CT is the absolute value of the electronic coupling between the CT and the ground state. To determine the donor–acceptor electronic coupling, a statistically representative ensemble of donor–acceptor configurations must be considered. Indeed, in disordered organic semiconductors, the coupling is highly sensitive to the relative donor–acceptor geometry. Within our framework, classical molecular dynamics (MD) simulations are employed to sample local donor–acceptor interfacial arrangements and generate a statistically representative ensemble of configurations.^79,80 The electronic coupling is then computed for each sampled configuration according to:

V_CT = |〈ϕ_h|F̂|ϕ_ℓ〉|, (34)

where ϕ_h and ϕ_ℓ are the unperturbed HOMO of isolated D and the unperturbed LUMO of isolated A, respectively, and F̂ is the Kohn–Sham-Fock operator of the D:A complex.^81–83 Finally, an effective coupling is obtained as a statistical average over the ensemble (see Computational details). This approach has been successfully applied in numerous theoretical studies of charge-transfer states in organic photovoltaic blends.^{13,14,84–86}

Although real bulk heterojunction films can exhibit complex mesoscale features, such as vertical phase segregation, variations in domain purity, and composition gradients that depend on processing conditions, these features mainly affect the amount and spatial distribution of donor–acceptor interfaces, rather than the electronic interactions at an individual donor–acceptor contact. Since charge carriers are assumed to have infinite mobility within the present framework (see above), mesoscale transport pathways do not explicitly enter the model. Accordingly, molecular dynamics simulations, used only to infer V_CT, are aimed at sampling local donor–acceptor interfacial configurations relevant for charge-transfer state formation rather than reproducing the full mesoscale morphology of bulk heterojunction films. Therefore, the evaluation of V_CT based on local interfacial configurations is sufficient for capturing the physics of the CT state addressed in this work.

Results and discussion

The method was tested on the binary cells experimentally investigated and listed in Table 1, which reports each D:A couple along with the experimental weight ratio, the measured and predicted open-circuit voltage, short-circuit current density and the observed power conversion efficiency. The molecular structures of the investigated donor and acceptor materials are provided in Fig. S1 and S2 in the SI.

Table 1 Donor:acceptor (D:A) pairs together with their mass ratio and observed power conversion efficiency (PCE). Observed (exp) and predicted (theor) open-circuit voltage (V_OC, V) and short-circuit current density (J_SC, mA cm⁻²). All data refer to T = 300 K

D:A	Weight ratio	V^(exp)OC	V^(theor)OC	J^(exp)SC	J^(theor)SC	PCE (%)
a Ref. 87.b Ref. 88.c Ref. 39. J^(exp.)SC refers to the cell with composition 1 : 19 ratio.d Ref. 16.e Ref. 89.f Ref. 90.
FG3:Y6^a	1 : 1.5	0.90	0.95	18.38	27.39	10.75
FG4:Y6^a	1 : 1.5	0.79	0.84	20.92	24.98	11.07
DTS:PCBM70^b	1 : 0.7	0.80	0.80	13.40	13.06	6.02
α-6T:C60^c	1 : 1.5	0.46	0.55	5.70	7.30	0.90
ZR1:Y6^d	1 : 0.5	0.86	0.89	24.34	27.44	14.27
ZR1:IDIC-4Cl^d	1 : 0.7	0.78	0.87	18.27	25.59	9.58
BTEC-2F:Y6^e	1 : 1	0.85	0.90	21.55	24.33	13.34
DRCN5T:PCBM70^f	1 : 0.8	0.92	0.90	15.66	12.60	9.80

---

The selected blends display V_OC values in the range 0.46–0.92 V and J_SC values between approximately 6 and 24 mA cm⁻², with PCE spanning from 0.9 to 14.3% (Table 1). These ranges cover a large portion of experimental results reported in the literature, reflecting the typical performance window observed for bulk heterojunction solar cells across various donor–acceptor compositions.

As shown in Table 1 (see also Fig. S2 in the SI), we considered both fullerene (C₆₀, PCBM70) and non fullerene acceptors (NFA) such as Y6 and IDIC-4Cl, which address key limitations of fullerene based systems, including low absorption coefficients, poor chemical stability, tendency to self aggregate into pure domains that reduce the donor acceptor interfacial area, and limited energy tunability.^16,91

The donor structures are shown in Fig. S1 of the SI. FG3 and FG4 were used for their complementary absorption spectra with Y6 as the acceptor.^31,87 Additional donors include ZR1, known for short π−π stacking that enhances hole transport,¹⁶ DRCN5T, an A–D–A donor with an oligothiophene core,⁹⁰ DTS, a D–A–D type low band-gap donor available commercially,⁸⁸ BTEC-2F, featuring fluorinated thiophene units,⁸⁹ and α-6T, composed of six thiophene rings, often used as a model compound able to capture the key electronic features of conjugated polymeric donors.⁹²

The structural characteristics of both donors and acceptors are reflected in the FCWD. The latter plays a pivotal role also because it permits to infer reliable estimates of the ΔE₀₀ adiabatic energies associated with the S₁ ← S₀ excitation of the cell components.

Fig. 2 compares the FCWDs computed via the GF approach with those experimentally extracted from the absorption spectra of the BTEC-2F donor and IDIC-4Cl acceptor recorded in thin films.^16,89 Theoretical profiles reproduce the experimental band shapes and widths without any additional broadening or fitting parameters.

Fig. 2 Predicted (blue) and observed (red) FCWDs for the S₁ ← S₀ absorption of BTEC-2F and IDIC-4Cl at T = 300 K (top and bottom, respectively).^16,89 Vertical dotted lines mark the adiabatic excitation energy (ΔE₀₀) values obtained by spectral alignment. The experimental FCWDs extend further to the right due to additional absorptions leading to higher energy excited states not included in the calculations.

In our methodology, the purely electronic gap ΔE is initially set to zero in eqn (42), and the computed FCWD is rigidly shifted along the energy axis to best match the experimental FCWD; this energy shift directly yields the ΔE₀₀ value.

Previous studies have shown that the ΔE₀₀ values obtained by aligning the computed FCWD with the experimental FCWD of either absorption or emission (even phosphorescence) provides more accurate ΔE₀₀ estimates than those derived solely from DFT or TDDFT computations (see also Table S1 in the SI and ref. 68 and 75). The ΔE₀₀ values obtained through this procedure, reported in Table 2 as E_D and E_A, are then used in eqn (18) to compute α for all donors and acceptors, except for C₆₀ and PCBM70, whose adiabatic excitation energies are taken from ref. 93, due to the well-known difficulties in resolving the S₁ ← S₀ absorption band of fullerenes.⁹⁴

Table 2 Excitation energies of the donor and acceptor (E_D, E_A, eV) obtained by aligning the experimental and computed FCWDs for S₁ ← S₀ absorption, except for C₆₀ and PCBM70, see main text. Charge-transfer energies (E_CT, eV) from voltammetric experiments, unless otherwise specified. Predicted dipole strengths for the donor (|μ^D_tr|²) and the acceptor (|μ^A_tr|²) as isolated species. Dipole strength for the CT state (|μ^CT_tr|²) predicted by using the diabatization method. All dipole strengths are given in Debye². Predicted electronic couplings (V_CT, meV) used in the calculation of non-radiative recombination rates

D:A	ED	EA	ECT	|μ^Dtr|^2	|μ^Atr|^2	|μ^CTtr|^2	VCT
a Ref. 87.b Ref. 88.c Ref. 93.d Ref. 95.e |μ^Atr|^2 obtained by the diabatization method.f Ref. 96.g The ionization potential of α-6T has been determined via photoelectron spectroscopy on 3D thin films,^97 whereas the oxidation free energy of C60 has been inferred from voltammetry measurements using ferrocenium/ferrocene (Fc^+/Fc) as internal reference and carried out on thin films of C60 adhered to Pt electrode in acetonitrile;^98 the absolute reduction potential of the Fc^+/Fc couple has been set at 4.98 V.^99h Ref. 16.i Ref. 89.j Ref. 90.
FG3:Y6	1.65^a	1.48^a	1.33^a	603.66	385.18	125.79	4.9
FG4:Y6	1.48^a	1.48^a	1.11^a	842.91	385.18	138.50	4.1
DTS:PCBM70	1.79^b	1.80^c	1.21^b^d	337.87	4.75^e	7.69	7.0
α-6T:C60	2.37^f	1.85^c	1.03^g	235.65	0.15^e	0.17	14.1
ZR1:Y6	1.95^h	1.48^a	1.26^h	506.86	385.18	22.17	6.6
ZR1:IDIC-4Cl	1.95^h	1.57^h	1.22^h	506.86	479.47	249.57	4.8
BTEC-2F:Y6	1.95^i	1.48^i	1.29^i	369.50	385.18	166.41	19.3
DRCN5T:PCBM70	1.79^j	1.80^c	1.31^d^j	425.94	32.9^e	58.43	10.5

---

Table 2 also presents, for each D:A pair, the dipole strength associated with the S₁ ← S₀ transition of all donor and acceptor species, as predicted by TDDFT calculations. For the weakly absorbing C₆₀ and PCBM70 systems, the dipole strength was obtained via the diabatization procedure outlined in the Computational details. The weak absorption of fullerene acceptors (Table 2) is also reflected in the low dipole strength of the CT states formed with the donor.

The CT energies in Table 2 were calculated from eqn (21) using absolute redox potentials of donors and acceptors measured by voltammetry in the same laboratory, whenever possible (see Table S1 in the SI). For α-6T, voltammetric data are unavailable because it tends to polymerize on commonly used electrodes.^100,101 Thus, in lieu of the absolute oxidation potential, its ionization potential was used, as obtained via photoelectron yield spectroscopy on 3D thin films.⁹⁷

Table 2 shows that CT energies range from approximately 1.0 to 1.3 eV and couplings between the ground state and CT state (V_CT) are of the order of tens of meV, as expected.²⁸

Predicted V_OC and J_SC are reported in Table 1 and plotted in Fig. 3 against measured values.

Fig. 3 Predicted vs. observed V_OC (top) and J_SC (bottom) at T = 300 K. Solid lines have null intercept and unitary slope.

Comparison of predicted and observed V_OC values (Table 1 and the top panel of Fig. 3) indicates generally good agreement. The mean absolute error (MAE) amounts to 0.05 V, corresponding to an average percent error of about 6%, with deviations within 0.05 V for six of the eight devices examined. The mean signed error (MSE) is +0.04 V, reflecting a slight overestimation.

The predicted values of J_SC (Table 1 and bottom panel of Fig. 3) yield an MAE of 3.89 mA cm⁻². The quantitative agreement for J_SC is overall satisfactory within the present level of approximation. At the same time, the model correctly discriminates between fullerene- and non-fullerene-based D:A blends through the predicted J_SC values, yielding systematically lower J_SC for fullerene-based cells than for non-fullerene-based ones, in line with their different optical absorptance (eqn (11)).

The MSE is 3.06 mA cm⁻², indicating an overall positive deviation. Overestimated J_SC are found for blends employing non-fullerene acceptors. This behavior is consistent with the assumptions adopted here, namely an infinite-mobility limit and the neglect of trap-assisted recombination, as fullerene-based bulk heterojunctions are typically found to exhibit higher mobilities than NFA-based cells.^102–104 Further evidence supporting this analysis is provided by comparing the J_SC values of the FG3:Y6 and FG4:Y6 systems (Table 1). For FG3:Y6, the computed J_SC exceeds the experimental value by approximately 9 mA cm⁻², whereas for FG4:Y6 the overestimation is reduced to only 4 mA cm⁻². Under the same assumptions of infinite carrier mobility and neglected trap-assisted recombination, the different magnitude of this overestimation can be rationalized by the higher hole mobility and weaker trap-assisted recombination in FG4:Y6 compared to FG3:Y6, as determined from experimental studies.⁸⁷

An overestimation of J_SC could, in principle, lead to an overestimation of V_OC,rad (eqn (10)), which in turn would suggest that the slight overestimation of V_OC originates from an overestimated radiative component. However, this may not be the case. For the α-6T:C₆₀, ZR1:Y6, and ZR1:IDIC-4Cl devices, experimental determinations of both V_OC and ΔV_OC,nr are available (Table 3). In these cases, the comparison between predicted and measured ΔV_OC,nr reveals that the slight overestimation of V_OC mainly arises from an underestimation of ΔV_OC,nr. Considering, for instance, the blend exhibiting the largest deviation in V_OC, namely α-6T:C₆₀ (0.09 V), the predicted V_OC is 0.55 V, while the observed value is 0.46 V. For this blend, the predicted nonradiative voltage loss ΔV_OC,nr amounts to 0.41 V, compared to the observed value of 0.46 V, thereby confirming an underestimation of ΔV_OC,nr by the model. The occurrence of the largest deviation for the α-6T:C₆₀ cell is somewhat expected, as the donor IP and acceptor EA entering E_CT (eqn (21)) are not fully consistent. Specifically, the IP of α-6T (4.84 eV) is derived from photoelectron spectroscopy on 3D molecular sheets,⁹⁷ whereas the EA of C₆₀ (3.81 eV) is inferred from the reduction free energy measured by voltammetry on films formed on platinum electrodes in acetonitrile.⁹⁸ A more internally consistent estimate could be obtained by employing the solid-state EA of C₆₀. Indeed, assuming for C₆₀ an EA of 3.90 eV, well compatible with the range of reported values 3.5–4.4 eV,^105,106 instead of 3.81 eV would yield an almost perfect match between predicted and experimental V_OC. This highlights the need for accurate determinations of E_CT for quantitative predictions.

Table 3 Predicted values of V_OC,rad, ΔV_OC,nr (V); k_rad and k_nr (s⁻¹) at T = 300 K

D:A	VOC,rad	ΔVOC,nr	krad	knr
a Observed ΔVOC,nr = 0.46 V.^39b Observed ΔVOC,nr = 0.24 V.^16c Observed ΔVOC,nr = 0.38 V.^16
FG3:Y6	1.11	0.16	5.84 × 10^8	3.61 × 10^10
FG4:Y6	1.00	0.16	4.45 × 10^8	5.68 × 10^10
DTS:PCBM70	1.06	0.26	3.74 × 10^7	1.91 × 10^11
α-6T:C60	0.96	0.41^a	4.23 × 10^5	7.56 × 10^11
ZR1:Y6	1.09	0.20^b	5.60 × 10^7	2.17 × 10^10
ZR1:IDIC-4Cl	0.99	0.12^c	6.00 × 10^8	1.44 × 10^10
BTEC-2F:Y6	1.07	0.17	4.90 × 10^8	6.42 × 10^10
DRCN5T:PCBM70	1.11	0.21	2.81 × 10^8	2.06 × 10^11

---

Table 2 shows that, for NFA-based cells, TDDFT calculations generally predict large dipole strengths for the charge-transfer state, comparable in magnitude to those of the individual components. These values appear to be overestimated, since their use leads to predicted maximum values of α_CT that occasionally exceed ∼10³ cm⁻¹, as shown in Fig. 4, whereas the measured absorption coefficients in the subgap region typically amount to only a few hundred cm⁻¹.⁶⁴ Therefore, to assess the impact of |μ^CT_tr|² on V_OC,rad, ΔV_OC,nr, and J_SC, we performed numerical tests in which the dipole strength was scaled by factors of 1/25 and 1/100. As shown in Table 4, upon reducing the oscillator strength, for NFA-based cells both J_SC and V_OC,rad remain nearly constant and very close to the values obtained using the unscaled dipole strengths (compare with Table 3). Thus, even if CT dipole strengths are somewhat overestimated, this does not significantly affect the predicted absorptance. On the opposite, ΔV_OC,nr increases by about 0.09 V when |μ^CT_tr|² is reduced by a factor of 25, and by about 0.12 V when it is lowered by a factor of 100. This provides, although not conclusive, a further indication that predicted ΔV_OC,nr are possibly slightly underestimated due to overestimated dipole strengths for emission. The present implementation employs two approximations. First, the same dipole strength is used for both absorption and emission, which, as anticipated, is not expected to be critical. Second, unlike V_CT, which is averaged over several donor–acceptor configurations, |μ^CT_tr|² was evaluated only for the minimum-energy D:A complex, due to the computational cost of TDDFT. This second approximation may be more relevant, because averaging over multiple configurations of the D:A complex could significantly alter the dipole strength. In conclusion, our analysis suggests a slight overestimation of the dipole strength for the CT state, which in turn leads to a moderate overestimation of k_rad and consequently to a slight underestimation of ΔV_OC,nr. However, it remains unclear whether this effect originates from the lack of statistical sampling of |μ^CT_tr|², or from overestimation by the adopted functional (CAM-B3LYP, see the Computational details section).

Fig. 4 Absorption coefficient of the ZR1 donor (α_D, blue), IDIC-4Cl acceptor (α_A, red), and the charge-transfer state (α_CT, green), along with their sum (α, black dashed line), as a function of the incident light wavenumber, for the ZR1:IDIC-4Cl device at T = 300 K. The spectral photon flux of the black body (ϕ_BB) is shown in purple. The inset highlights the spectral region where α_CT, despite its comparatively weak intensity relative to the other components, contributes to V_OC.

Table 4 Open-circuit voltage components (V_OC,rad and ΔV_OC,nr, V) and short-circuit current density (J_SC, mA cm⁻²) calculated at T = 300 K after dividing the squared transition dipole moment |μ^CT_tr|² by the factors f = 25 and f = 100

D:A	VOC,rad		ΔVOC,nr		JSC
	f = 25	f = 100	f = 25	f = 100	f = 25	f = 100
FG3:Y6	1.11	1.11	0.24	0.27	27.39	27.39
FG4:Y6	1.07	1.08	0.25	0.29	24.90	24.89
ZR1:Y6	1.10	1.10	0.28	0.32	27.43	27.43
ZRI:IDIC-4Cl	1.08	1.12	0.20	0.24	25.28	25.27
BTEC-2F:Y6	1.09	1.09	0.26	0.30	24.31	24.30

---

The impact of direct formation of the CT state on the open-circuit voltage can be argued by the inspection of Fig. 4, which displays the predicted total absorption coefficient of the ZR1:IDIC-4Cl blend as a function of the incident light wavenumber, together with the individual contributions α_A, α_D, and α_CT. In the range of validity of the thin-film limit (eqn (16)), the absorptance, i.e., the external quantum efficiency (eqn (14)), is proportional to α(E). Since the product of the absorptance and ϕ_BB enters the expression of J_0,rad (eqn (12)), the larger the overlap between α(E) and ϕ_BB, the larger J_0,rad. Therefore, if the absorption leading to charge-transfer states occurs at low energy, where ϕ_BB is large, J_0,rad is large as well. This, in turn, leads to a decrease of V_OC,rad according to eqn (10). Thus, although α_CT is much weaker than α_D and α_A, its contribution may still adversely affect the V_OC of the device, if it falls at low energy. That detrimental effect was indeed observed both in BHJ and perovskite-based solar cells.^48,65 Conversely, if the absorption onset occurs at higher energies, the exponential decrease of ϕ_BB (purple line in Fig. 4) keeps the product α(E)ϕ_BB small. As a result, the overlap remains limited, and J_0,rad, and hence V_OC,rad, are only weakly affected. In addition, ΔV_OC,nr varies only slightly when α_CT is included in the prediction of α(E), see eqn (29) and (31), hence the total V_OC decreases by roughly the same amount as the reduction in V_OC,rad upon considering α_CT. As a numerical test, omitting α_CT in the ZR1:IDIC-4Cl calculations shifts V_OC,rad from 0.99 to 1.21 V, while ΔV_OC,nr increases only marginally, from 0.12 to 0.14 V.

Present theoretical predictions for V_OC not only show good agreement with experimental data, but also reflect the correct behavior of organic solar cells and the nature of its costituents. Notably, for the absorption of fullerene-based acceptors, whether the transition dipole moment |μ^A_tr| is set to zero (as is the case for the lowest absorption of I_h-symmetric C₆₀) or evaluated via the diabatization procedure, the resulting variations in the predicted open-circuit voltage, V_OC, are below approximately 0.02 V, while changes in the short-circuit current density, J_SC, remain under ≈0.5 mA cm⁻².

Moreover, V_OC,rad is correctly calculated to increase linearly with E_CT.³⁹ This is evident in the upper panel of Fig. 5, where computed V_OC,rad is plotted against E_CT, yielding a clear linear trend with R² = 0.79. The bottom panel of the same figure shows the total V_OC as a function of the energy offset ΔE_LE–CT, i.e., the driving force, herein defined as the difference between the lower of the donor and acceptor excitation energy and E_CT. As predicted by our model, and consistent with simple energetic considerations, a larger offset inevitably leads to additional losses and thus a lower V_OC.⁴⁰

Fig. 5 Top: Predicted V_OC,rad as a function of E_CT. The solid line is obtained from a linear fit. Bottom: predicted V_OC as a function of the energy offset, ΔE_LE–CT, defined as the difference between the lower of the donor and acceptor excitation energy and E_CT. Red diamonds denote fullerene based cells; blue circles NFA-based cells. The temperature is set to 300 K.

Table 3 shows that for cells with fullerene acceptors, k_nr, i.e., the non-radiative decay rates of the CT state, are predicted to be generally larger than those for non-fullerene acceptors, as indeed observed experimentally.¹⁰⁷ This increase in non-radiative recombination rates correlates with higher values of predicted ΔV_OC,nr in fullerene-based devices compared to non-fullerene systems, in agreement with experimental observations. Indeed, non fullerene acceptors have been applied in organic solar cells also because they ensure lower voltage losses, caused by lower non radiative recombination rates, when compared to fullerene based systems.^107–109 The top panel of Fig. 6 depicts the non-radiative voltage loss as a function of E_CT. Fullerene-based acceptor cells (red diamonds), though only three data points are available, exhibit a clear decreasing linear trend. This pattern was indeed observed for organic solar cells with fullerene acceptors³⁹ and is consistent with the energy gap law.^110,111 In contrast, blends based on non-fullerene acceptors (blue circles) exhibit no correlation between ΔV_OC,nr and E_CT, as indeed observed.²⁸ This lack of correlation is supported by the low R² = 0.33 value obtained when considering all data points. The same conclusions emerge from the bottom panel of Fig. 6, where ΔV_OC,nr is plotted as a function of ΔE_LE–CT. Blends based on non-fullerene acceptors (blue circles) again show no clear trend. In contrast, fullerene-based cells (red diamonds) exhibit a clear linear decrease. This trend, too, is consistent with the energy gap law, which has been demonstrated to apply only to fullerene-based cells, provided that ΔE_LE–CT ≥ 0.2 eV,²⁸ a value well below those in our data set, where ΔE_LE–CT systematically exceeds 0.4 eV.

Fig. 6 Predicted ΔV_OC,nr as a function of E_CT (top) and ΔE_LE–CT (bottom) at T = 300 K. Red diamonds denote fullerene based cells; blue circles NFA-based cells. The solid line is obtained from a linear fit.

Since our treatment of the Franck–Condon weighted density of states intrinsically accounts for thermal effects, it also allows the temperature dependence of the open-circuit voltage and short-circuit current density to be assessed, which is essential for characterizing device performance under realistic outdoor operating conditions. The assessment covers three of the investigated systems, namely α-6T:C₆₀, DTS:PCBM70, and BTEC-2F:Y6. The results are shown in Fig. 7, which reports V_OC and J_SC as a function of temperature in the range 150–400 K, the latter taken as an upper bound for model validation rather than steady-state device operation, and are detailed in Table S2 of the SI, where the corresponding radiative and nonradiative decay rates, k_rad and k_nr, and the radiative and nonradiative contributions to V_OC at the temperatures considered are also reported.

Fig. 7 Open-circuit voltage (top) and short-circuit current density (bottom) plotted as a function of the temperature for α-6T:C₆₀ (left), DTS:PCBM70 (center), and BTEC-2F:Y6 (right). The red circles denote the computed values, while the dashed grey lines are obtained through linear fitting.

Inspection of the top panels of Fig. 7 shows that the open-circuit voltage is predicted to decrease linearly with increasing temperature. That linear trend has been observed in similar systems and has been attributed to two main effects: (i) a broadening of the absorption onset with increasing temperature,⁵¹ and (ii) an enhancement of nonradiative recombination losses.^{51,52,112,113} Within our framework, the temperature dependence of the FCWD results in a broader absorptance which arises because a larger number of vibrational modes becomes thermally activated as temperature increases. This effect, via eqn (12), leads to an increased overlap of absorptance with the black-body photon flux. Since the black-body spectrum also broadens with temperature, this causes the radiative saturation current J_0,rad to increase more rapidly than J_SC. As a result, the reduction of the logarithmic term (eqn (10)) outweighs the explicit temperature prefactor k_BT/q, leading to a decrease of V_OC,rad with increasing temperature. The second effect is accounted for by the temperature dependence of charge-transfer state decay: both k_rad and k_nr are predicted to increase with temperature, with k_nr exhibiting a steeper temperature dependence. This, in turn, causes an increase in ΔV_OC,nr (eqn (29)), thus further reducing V_OC upon increasing temperature, with the same considerations as before for the temperature prefactor.

Turning to a quantitative comparison, the temperature derivatives of the open-circuit voltage predicted by our model for the α-6T:C₆₀, DTS:PCBM70, and BTEC-2F:Y6 devices are −1.7, −1.2, and −1.6 mV K⁻¹, respectively. Notably, for the planar heterojunction α-6T:C₆₀ device, a temperature derivative of −1.9 mV K⁻¹ can be inferred from the experimental data reported in ref. 112. Although temperature derivatives are not available for the other systems investigated here, values of −1.3 mV K⁻¹ for the DIBSQ:PCBM70 BHJ and −2.9 mV K⁻¹ for the planar heterojunction DBP:C₇₀ have been reported from experimental measurements in ref. 51 and 52, respectively. Overall, these comparisons indicate that the present model not only captures the qualitative trend of V_OC with temperature, but also provides quantitatively reliable estimates of its temperature dependence.

The bottom panels of Fig. 7 show that our model predicts a linear increase of J_SC with temperature. Available experimental studies support this trend, reporting an increase of J_SC with temperature,^51,52,113 which is found to be approximately linear over a broad temperature range around 300 K.^51,52 This behavior is generally attributed to two main effects, both enhanced by increasing temperature: (i) the broadening of the absorptance A(E), which increases its overlap with the solar spectrum (eqn (11)), and (ii) the increase of the charge-carrier mobility.

Since our model assumes infinite carrier mobility, the predicted increase of J_SC originates solely from the broadening of A(E) and the resulting enhanced overlap with ϕ_sun. This effect is shown in Fig. S3 of the SI. In line with this interpretation, for the planar heterojunction DBP:C₇₀, the short-circuit current density was observed to increase linearly with temperature at a rate of 6.0 × 10⁻³ mA cm⁻² K⁻¹ between 298 and 398 K, and this behavior was attributed primarily to absorptance broadening.⁵² In other systems, such as the DIBSQ:PCBM70 BHJ, a temperature derivative of 1.9 × 10⁻² mA cm⁻² K⁻¹ has been measured, where the dominant contribution was attributed to the increase in carrier mobility.⁵¹ Despite the simplifying assumptions of our model, we obtain slopes of 7.4 × 10⁻³, 9.9 × 10⁻³, and 2.1 × 10⁻² mA cm⁻² K⁻¹ for α-6T:C₆₀, DTS:PCBM70, and BTEC-2F:Y6, respectively. These values are consistent with experimentally reported temperature derivatives and demonstrate that the model also provides realistic estimates of the temperature dependence of J_SC.

The present approach employs four parameters: the refractive index (default n = 1.5), the film thickness (default d = 100 nm), the blend density (default ρ = 1.2 kg dm⁻³), and the ratio of the number density of the active CT complex and the minimum of the number densities of the donor and the acceptor (default r = 10⁻²). In general, the model exhibits low sensitivity to these parameters, and the variations (at T = 300 K) resulting from their modification follow the expected trends.

When increasing the refractive index from n = 1 to n = 2, V_OC (Table S3 in the SI) is predicted to change by less than 0.01 V; small variations are also predicted for J_SC, which increases with n, consistent with enhanced absorption (see eqn (18)).⁴⁵ However, this increase, about 2 mA cm⁻², is smaller than the mean absolute error between calculated and observed J_SC values at n = 1.5. Furthermore, predicted p_e (Table S4 in the SI) is in reasonable agreement with eqn (31).

Our predicted V_OC exhibits minimal sensitivity to active-layer thickness, decreasing by roughly 0.03 V, which is smaller than the MAE of V_OC at d = 100 nm, when the thickness increases from 50 to 200 nm (Table S5 in the SI). J_SC increases on average by 5.7 mA cm⁻² over the same range, in close agreement with experimental observations.^114–116

The same trend for V_OC, namely an average decrease of about 0.04 V (see Table S6 in the SI), is predicted by our model when increasing the active-layer density from 0.8 to the limiting value of 1.6 kg dm⁻³. In this range J_SC is predicted to increase, on average by 3 mA cm⁻², as a result of the enhanced absorption associated with the higher density, supporting the same considerations made for the thickness dependence.

Regarding the number density ratio r of the active CT complex (Table S7 in the SI), when varying r from 10⁻³ to 10⁻¹, the predicted J_SC changes negligibly, with an average increase of only 0.9 mA cm⁻². Comparison of data in Table 3 and Table S7 in the SI shows that V_OC decreases in two nearly equal steps, each averaging about 0.045 V, as r increases from 10⁻³ to the default value 10⁻², and then from 10⁻² to 10⁻¹, thus exhibiting a logarithmic trend. In particular, upon raising r from the default (10⁻²) to the unrealistically large 10⁻¹, the average variation of V_OC matches its mean absolute error evaluated at r = 10⁻², thus indicating, on the whole, a moderate sensitivity of V_OC to r. Because ΔV_OC,nr and J_SC remain essentially unchanged with respect to r, the effect of r concerns only V_OC,rad, which depends logarithmically on the inverse of the dark-saturation current J_0,rad (eqn (10)). As a result, the increase in r leads to an increase in J_0,rad, thus lowering V_OC,rad. This clearly reflects the effect illustrated in Fig. 4, which shows that the absorption due to the CT complex – amplified by increasing r – adversely affects V_OC.

Table S8 in the SI reports V_OC and J_SC computed with the excitation energies of D and A, as well as the CT energies obtained from DFT/TDDFT calculations using the B3LYP functional. When compared with experimental data, it emerges that for J_SC, both the MAE and the MSE are approximately doubled relative to predictions that employ experimental excitation and CT energies. The deviation from the measured values is even more pronounced for V_OC, where the use of calculated energies yields an MAE of 0.15 V, compared to 0.05 V when experimental energies are adopted. The V_OC values obtained from calculated energies are generally overestimated (MSE = 0.09 V) because V_OC,rad are systematically overpredicted, while ΔV_OC,nr are underestimated. A striking example is provided by α-6T:C₆₀, which corresponds to the system with the largest deviation relative to experiment. For this case, the model based on DFT-derived energies predicts a V_OC overestimated by 0.53 V and a ΔV_OC,nr underestimated by 0.09 V.

Inspection of Tables S9 and S10 in the SI reveals that predicted excitation energies are, on average, underestimated by 0.27 eV, whereas CT energies are overestimated by 0.31 eV. This discrepancy mainly originates from a marked underestimation (−0.31 eV) of the electron affinity of the acceptors. These trends clearly highlight the well-known limitations of DFT/TDDFT methods in accurately predicting the energies of excited states and anions,^117–121 and they also imply that the predicted photovoltaic cell characteristics are severely distorted by such inaccuracies. Specifically, the energy gap between the donor absorption and the CT transition is strongly reduced, leading to an erroneous prediction of CT absorption at higher frequencies which causes an artificial increase in V_OC,rad. Furthermore, higher E_CT values generally correspond to regions where the FCWD for non-radiative emission is less populated. Consequently, the use of DFT energies results in underestimated k_nr and, therefore, in underestimated ΔV_OC,nr values. As an example, for the BTEC-4Cl:Y6 cell, when replacing the experimental E_CT with the DFT-predicted value, the FCWD evaluated at 300 K (eqn (33)) decreases from 2.24 × 10⁻⁶ to 8.75 × 10⁻⁸ cm, lowering k_nr from 6.4 × 10¹⁰ to 2.5 × 10⁹ s⁻¹. As a result, ΔV_OC,nr decreases from 0.17 V (Table 3) to 0.08 V (Table S8, SI).

Conclusions

We have developed a novel computational framework to investigate the factors governing the open-circuit voltage in bulk heterojunction solar cells, and to predict, with reasonable accuracy, both the open-circuit voltage and the short-circuit current density of small-molecule-based organic devices employing either fullerene or non-fullerene acceptors. The core of the method is the calculation of the external quantum efficiency of the cell and the rates of the elementary processes by a rigorous quantum mechanical approach, based on time dependent perturbation theory in which the Franck–Condon weighted density of states is computed from first principles, while MD simulations are used to sample donor–acceptor configurations in order to obtain averaged electronic couplings. The approach requires only a few experimental inputs, namely the redox potentials and absorption spectra of the neat donor and acceptor materials, and makes use of a small number of adjustable parameters whose values have only a minor effect.

Although the model is highly simplified, as it neglects many complexities present in real devices, such as trap-assisted and surface recombination, electrode effects and energetic disorder, it shows encouraging agreement with experimental data, demonstrating its capability to quantitatively predict both V_OC and J_SC. Its main strengths are: (1) low sensitivity to variations in adjustable parameters, (2) the ability to interpret device performance by attributing voltage losses to specific mechanisms, such as non-radiative recombination kinetics, subgap effects and the broadening of the absorption edge, nearly unavoidable phenomena in organic photovoltaics, and (3) the ability to quantitatively reproduce the temperature dependence of V_OC and J_SC, which is a key element for a meaningful assessment of the practical performance of bulk heterojunction solar cells under the variable operating temperatures encountered in outdoor conditions.

A primary limitation of the method as a purely predictive tool is its dependence on accurate estimates of the optical and fundamental gaps of both donor and acceptor. Experimental determination of these quantities remains at the present stage unavoidable.

Furthermore, in its present form, the method is not applicable to conjugated polymeric systems, which require dedicated theories of optical transitions and exciton transfer.^122–124

Of note, the results are remarkably good given the simplicity of the model, the more so as further refinement is possible e.g., by improving the prediction of dipole transition moment of the CT state, introducing a more realistic trapping scheme that includes Lambertian scattering at the front surface,¹²⁵ and by relaxing the assumption of virtually infinite charge-carrier mobility.⁶¹ In this direction, we plan to extend the current approach by incorporating a kinetic Monte Carlo framework that explicitly accounts for charge transport and mobility effects.^14,126–128

Computational details

MD simulations

Classical molecular dynamics simulations were performed on a cubic box comprising a number of molecules per species selected to reflect the experimental blend composition. The molecules have been randomly arranged to mimic a BHJ interface. Computations were carried out by using the GROMACS 2020.5 software,¹²⁹ following the procedure reported in ref. 13. The OPLS force field was adopted throughout.¹³⁰ Periodic boundary conditions were introduced and long range electrostatic effects were taken into account through the particle mesh Ewald algorithm (PME).¹³¹ An integration time step of 2 fs was set up, by imposing constraints on the bonds involving H atoms through the linear constraint solver (LINCS) algorithm.¹³² A modified Berendsen thermostat was adopted to control temperature and a Parrinello–Rahman barostat to control pressure.^133,134 The computational protocol consisted of an initial steepest descent minimization, followed by a NVT equilibration at 300 K over 1 ns. The system was then subjected to 5 simulated annealing (SA) cycles over 10 ns, between 300 K and 600 K, adopting the NPT ensemble: heating and cooling steps were run over 500 ps intervals, the target temperature was maintained for 9.5 ns. Finally, a production run of 100 ns in the NPT ensemble was carried out.

Calculation of charge transfer couplings

In order to compute the electronic couplings V_CT, for each cell, at least 100 D:A pair geometries were extracted from the last 40 ns of the simulation. These configurations are statistically independent, being separated by at least 500 ps. Additional convergence tests were performed for two representative systems, ZR1:Y6 and ZR1:IDIC-4Cl, for which V_CT was averaged over a much larger number of D/A pairs. For these two blends, the mean values of V_CT obtained by averaging over an increasing number of D/A pairs are reported in Fig. S4 of the SI. In both cases, about 100 configurations were found to be sufficient to yield stable mean values and narrow distributions of the electronic couplings. Therefore, at least 100 configurations were used for the remaining blends. Following ref. 13, only configurations in which the D and A molecules had a long molecular axis at least 90% as long as in their equilibrium geometry, and with at least two atoms (one from D, one from A) closer than 5 Å, were retained. The required Kohn–Sham matrix elements (eqn (34)) were then computed at the density functional level of theory (DFT) in conjunction with the 3-21G* basis set, using the B3LYP functional.^31,135–138 That choice has been made to reduce the computational cost, given the weak dependence of the couplings on the basis set.⁸²

DFT and TDDFT computations

DFT and TDDFT calculations on isolated donor and acceptor molecules in their neutral ground (S₀) and first excited (S₁) states were performed using the B3LYP functional. For all donor and acceptor species considered here, the lowest-energy transition S₁ ← S₀ is an intense, nearly pure HOMO–LUMO excitation. Both π frontier orbitals are mostly delocalized across the entire conjugated backbone (see Fig. S5–S13 in the SI). This, along with the small change in permanent electric dipole moment upon excitation (Table S11 in the SI), indicates a valence-type transition which, as is well known, is reliably described by global hybrid functionals such as B3LYP.^139,140 Furthermore, it has been shown that TD-B3LYP predicted oscillator strengths for the HOMO–LUMO absorption of non-fullerene acceptors agree reasonably well with experimentally determined maximum molar extinction coefficients.⁴⁵ The same functional was also employed to predict the equilibrium geometry and Hessian matrix of cationic (D⁺) and anionic (A⁻) species required for FCWD calculations for the CT complex.

CAM-B3LYP was employed in TDDFT calculations of the transition dipole moments for the absorption leading to CT and used in the diabatization scheme. This change of functional was necessary since global hybrid functionals with low percentage of Hartree–Fock exchange, such as B3LYP, are known to produce spurious, dark low-energy transitions for charge-transfer.¹⁴¹ Recent studies indicate that this issue is particularly severe for intermolecular CT, which is even more difficult to model than intramolecular CT. For the former, only range-separated double-hybrid (RS-DH) functionals based on the algebraic diagrammatic construction method at second order, ADC(2), can provide excitation energies and oscillator strengths of good quality, comparable to high-level coupled cluster reference data.^142,143 RS-DH functionals are, however, computationally too demanding for the molecular sizes investigated in this work. Therefore, we opted for CAM-B3LYP, a long-range corrected functional, which is known to mitigate the self-interaction error underlying failures in TDDFT computations for CT transitions.^140,144,145 Furthermore, the CAM-B3LYP functional has been extensively used in literature, including in recent studies,^45,146–148 to characterize excited-state properties of donor–acceptor systems typically used in organic semiconductors and bulk heterojunctions.

The unrestricted formalism was used for species with unpaired electrons. Dielectric effects were included in DFT and TDDFT calculations via the polarizable continuum model (PCM),¹²¹ with the static dielectric constant set to 4.0.¹⁰⁵ All DFT and TDDFT runs were carried with the Gaussian package, by using the 6-31+G(d,p) basis set.¹⁴⁹ The D3 scheme proposed by Grimme to include dispersion energy was used throughout.¹⁵⁰ Long lateral alkyl chains were included in MD simulations, but they were replaced by methyl groups in DFT and TDDFT computations.

Diabatization

To compute the transition dipole moments for the absorption leading to the CT state (Fig. 1a), we relied on the most stable D:A complex predicted by molecular dynamics simulations for each investigated cell. For all these complexes, we performed TDDFT calculations for the first three adiabatic (a) singlet states, which consistently exhibited non-negligible transition dipole moments in all tested cases.

Three diabatic (d) states, representing a local excitation on the donor, a local excitation on the acceptor, and a CT state, respectively, were obtained from the adiabatic states using the diabatization scheme defined in ref. 151. Given a dimer system with n adiabatic singlet excited states {ψ^(a)₁, ψ^(a)₂,…,ψ^(a)_n}, one can define an orthogonal transformation C relating adiabatic states to diabatic states:

As in many diabatization schemes, we seek a unitary transformation C that generates diabatic states close to reference states with well-defined character (e.g., Frenkel excitations localized on a single molecule). For most excited states, atomic transition charges allow monitoring the character of the excitation.¹⁵² We write:

q^(d,ref) = q^(a)C^T. (35)

Here, q^(a) is the row vector containing the atomic transition charges of the adiabatic state, while q^(d,ref) holds the reference transition charges constructed as follows:

• Transition charges for the first singlet excited state of the isolated donor (acceptor), setting to zero those on the acceptor (donor), to define the diabatic state localized on the donor (acceptor);

• Transition charges for the first singlet excited state of a positively charged donor and negatively charged acceptor, to define the CT state.

This distinguishes localized and CT-type excitations.

As discussed in ref. 151, a suitable transformation matrix C is defined by minimizing the Frobenius norm:

This corresponds to the orthogonal Procrustes problem,¹⁵³ with unique solution:

C^T = UV^T

where orthogonal U and V matrices are obtained from the singular value decomposition of

M = q^(a)Tq^(d,ref) = UΣV^T

with diagonal Σ. Once C is known, the transition dipole moment associated with the diabatic states can be computed. Simple algebra shows that, for each Cartesian component α = x, y, z, the following relation holds: μ^(d)_tr,α = μ^(a)_tr,αC^T,where μ^(d)_tr,α and μ^(a)_tr,α are row vectors containing the α component of the transition dipole moments of the diabatic and adiabatic states, respectively.

Prediction of the Franck–Condon weighted density of states

The kinetic rates and spectral line shapes appearing in the expression for α(ħω) are computed using the generating-function approach.⁷² The derivation and the implementation of the GF formalism were presented in previous work;^49,73 here, we provide only a brief overview of the main results.

According to Fermi's golden rule, the rate for an electronic transition from the initial state A to the final state B is given by:

where V_AB represents the proper electronic coupling between states A and B, as in eqn (32) and (33). Ω is the FCWD, and ΔE = E_B − E_A refers to the difference in “purely” electronic energies, with:Here, the plus (minus) sign in eqn (38) indicates emission (absorption).⁶⁶ By applying the Born–Oppenheimer approximation, for X = A, B, we have:

|X, x〉 = |X〉 ⊗ |x〉, (39)

E_X,x = E_X + E_x (40)

where |X〉 is the pure electronic state with energy E_X, and |x〉 represents one of the vibrational states of state X, whose energy is E_x. Assuming that V_AB is independent of nuclear coordinates (Condon approximation), the Franck–Condon weighted density of states (eqn (37)) can be rewritten as:⁷¹where β = 1/(k_BT), Z_A is the vibrational partition function of state A, and 〈b|a〉 are Franck–Condon integrals. Following Lax and Kubo, the FCWD can be expressed in terms of a generating function f(τ,T) as:^71,72where:in which is the vibrational Hamiltonian operator for the initial(final) state. Assuming harmonic vibrations and using mass-weighted vibrational coordinates, Kubo and Toyozawa were able to express the correlation function f(τ,T) in a closed form.⁷² More recently, our group extended and implemented the GF approach by employing mass-weighted normal coordinates in conjunction with Duschinsky's transformation.^{49,73,74,154,155} In order to compute f(τ,T), the equilibrium geometry and the Hessian matrix for the system in both states A and B are required; then the evaluation of the FCWD (eqn (42)) can be recast into a discrete Fourier transform problem.⁴⁹

All FCWD calculations were carried out using the freely available MolFC program.¹⁵⁵

The internal representation of normal coordinates was systematically adopted because it mitigates the occurrence of spurious couplings in the FCWD, which arise from unphysical shifts in stretching coordinates induced by the displacement of angular coordinates during the geometry change associated with the electronic transition. This issue is particularly pronounced when rectilinear coordinates, such as the Cartesian ones, are employed.¹⁵⁶

Owing to the intrinsic structural flexibility of the molecular components of BHJ cells, low-frequency vibrational modes are typically characterized by pronounced anharmonicity which can result in spurious features in the FCWD. While the most rigorous approach to remove such features would involve an explicit treatment of anharmonic effects, this entails a significant computational cost.^156,157 Hopefully, a similar goal can be achieved by excluding low-frequency modes from FCWD computations. This procedure, totally empirical in nature, has been demonstrated effective for comparable cases and is computationally much less demanding.⁷⁵ In this work, exclusion of all modes with harmonic wavenumbers below 100 cm⁻¹ proved sufficient for reliable results.

Recombination rates (eqn (33)) were computed by averaging the FCWD over a narrow interval (±100 cm⁻¹) around zero, a practical way to account for thermal disorder.³¹

Author contributions

The manuscript was written through the contributions of all the authors. All the authors have given approval to the final version of the manuscript.

Conflicts of interest

The authors declare no competing financial interest.

Data availability

The data supporting this article have been included as part of the supplementary information (SI). Supplementary information: molecular structures of donors (Fig. S1) and acceptors (Fig. S2). Absorptance weighted by solar irradiance as a function of wavenumber (Fig. S3). Averaged electronic couplings vs sampled configurations (Fig. S4). HOMO and LUMO isosurface contour plots for all molecules (Fig. S5–S13). Observed and predicted adiabatic absorption energies, donor ionization potentials, and acceptor electron affinities (Table S1). Predicted open-circuit voltage, short-circuit current density, and recombination rates at different temperatures (Table S2). Emission probability computed for different values of the refractive index (Table S4). Predicted open-circuit voltage and short-circuit current density as functions of refractive index (Table S3), device thickness (Table S5), device density (Table S6), and CT complex density (Table S7). Predicted open-circuit voltage, non-radiative losses, and short-circuit current density from PCM-TD-B3LYP-D3/6-31+G(d,p) excitation and CT energies (Table S8). Observed and predicted CT energies of D/A pairs (Table S9). MAE and MSE of excitation energy, CT energy, IP, and EA predicted by PCM-TD-B3LYP-D3/6-31+G(d,p) (Table S10). Predicted difference in permanent electric dipole moments between excited and ground states (Table S11). See DOI: https://doi.org/10.1039/d6qm00090h.

Acknowledgements

A. C., A. L. and A. P. gratefully acknowledge the support funding from the Italian Ministry of University and Research (MUR, PRIN grant 2022KHEZTC, 2022XSC9P5, and PRIN PNRR grant P2022WXPMB, respectively). T. P. acknowledges the PNRR funding from the Italian Ministry of University and Research (MUR, DM351− transizione digitale) for his PhD.

Notes and references

1. G. J. Hedley, A. Ruseckas and I. D. W. Samuel, Light Harvesting for Organic Photovoltaics, Chem. Rev., 2017, 117, 796–837.
2. J. Sun, X. Ma, Z. Zhang, J. Yu, J. Zhou, X. Yin, L. Yang, R. Geng, R. Zhu, F. Zhang and W. Tang, Dithieno[3,2-b:2′,3′-d]pyrrol Fused Nonfullerene Acceptors Enabling Over 13% Efficiency for Organic Solar Cells, Adv. Mater., 2018, 30, 1707150.
3. B. Walker, C. Kim and T.-Q. Nguyen, Small Molecule Solution-Processed Bulk Heterojunction Solar Cells, Chem. Mater., 2011, 23, 470–482.
4. Y. Sun, S. J. Edge, M. J. Ford, A. J. Heeger and G. C. Bazan, Solution-Processed Small-Molecule Solar Cells with 6.7% Efficiency, Nat. Mater., 2012, 11, 44–48.
5. S. E. Root, S. Savagatrup, A. D. Printz, D. Rodriquez and D. J. Lipomi, Mechanical Properties of Organic Semiconductors for Stretchable, Highly Flexible, and Mechanically Robust Electronics, Chem. Rev., 2017, 117, 6467–6499.
6. A. Landi, A. Peluso and A. Troisi, Quantitative Prediction of the Electro-Mechanical Response in Organic Crystals, Adv. Mater., 2021, 33, 2008049.
7. H. Zhang, Y. Li, X. Zhang, Y. Zhang and H. Zhou, Role of interface properties in organic solar cells: from substrate engineering to bulk-heterojunction interfacial morphology, Mater. Chem. Front., 2020, 4, 2863–2880.
8. B. D. Paulsen, K. Tybrandt, E. Stavrinidou and J. Rivnay, Organic Mixed Ionic-Electronic Conductors, Nat. Mater., 2020, 19, 13–26.
9. A. Landi, M. Reisjalali, J. D. Elliott, M. Matta, P. Carbone and A. Troisi, Simulation of polymeric mixed ionic and electronic conductors with a combined classical and quantum mechanical model, J. Mater. Chem. C, 2023, 11, 8062–8073.
10. G. P. Kini, S. J. Jeon and D. K. Moon, Latest Progress on Photoabsorbent Materials for Multifunctional Semitransparent Organic Solar Cells, Adv. Funct. Mater., 2021, 31, 2007931.
11. A. Köhler and H. Bässler, Electronic Processes in Organic Semiconductors, Wiley-VCH Verlag, Germany, 2015.
12. N. Yeh and P. Yeh, Organic Solar Cells: Their Developments and Potentials, Renewable Sustainable Energy Rev., 2013, 421–431.
13. A. Landi and D. Padula, Multiple Charge Separation Pathways in New-Generation Non-Fullerene Acceptors: a Computational Study, J. Mater. Chem. A, 2021, 9, 24849–24856.
14. A. Landi, A. Landi, A. Velardo and A. Peluso, Efficient Charge Dissociation of Triplet Excitons in Bulk Heterojunction Solar Cells, ACS Appl. Energy Mater., 2022, 5, 10815–10824.
15. J. Jin, Q. Wang, K. Ma, W. Shen, L. A. Belfiore, X. Bao and J. Tang, Recent Developments of Polymer Solar Cells with Photovoltaic Performance over 17%, Adv. Funct. Mater., 2023, 33, 2213324.
16. R. Zhou, Z. Jiang, C. Yang, J. Yu, J. Feng, M. A. Adil, D. Deng, W. Zou, J. Zhang, K. Lu, W. Ma, F. Gao and Z. Wei, All-Small-Molecule Organic Solar Cells with over 14% Efficiency by Optimizing Hierarchical Morphologies, Nat. Commun., 2019, 10, 5393.
17. T. Zhang, C. An, P. Bi, K. Xian, Z. Chen, J. Wang, Y. Xu, J. Dai, L. Ma, G. Wang, X. Hao, L. Ye, S. Zhang and J. Hou, A Highly Crystalline Donor Enables Over 17% Efficiency for Small-Molecule Organic Solar Cells, Energy Environ. Sci., 2024, 17, 3927–3936.
18. Y. Gao, L.-Y. Xu, X. Chen, B. Xiao, W. Gao, J. Xia, R. Sun and J. Min, Highly Efficient All-Small-Molecule Organic Solar Cells with Excellent Operational Stability and Blend-Thickness Tolerance, Energy Environ. Sci., 2025, 18, 7302–7312.
19. X. Cheng, Z. Liang, S. Liang, X. Zhang, J. Xu, Y. Xu, W. Ni, M. Li and Y. Geng, “Twisted” Small Molecule Donors with Enhanced Intermolecular Interactions in the Condensed Phase towards Efficient and Thick-Film All-Small-Molecule Organic Solar Cells, J. Mater. Chem. A, 2023, 11, 13984–13993.
20. C. Wang, T. Chen, S. Li, Y. Shen, J. Yu, A. Wupur, Y. Luo, M. Wang, X. Ye, J. Wu, M. Shi and H. Chen, Asymmetric Liquid Crystalline Donors with Two Different End Groups Enable Efficient All-Small-Molecule Organic Solar Cells, J. Mater. Chem. A, 2024, 12, 31163–31172.
21. H. Gao, Y. Sun, L. Meng, C. Han, X. Wan and Y. Chen, Recent Progress in All-Small-Molecule Organic Solar Cells, Small, 2022, 19, 2205594.
22. G. Han, Y. Zhang, W. Zheng and Y. Yi, Electron Transport in Organic Photovoltaic Acceptor Materials: Improving the Carrier Mobilities by Intramolecular and Intermolecular Modulations, J. Phys. Chem. Lett., 2023, 14, 4497–4503.
23. M. Riede, D. Spoltore and K. Leo, Organic Solar Cells-The Path to Commercial Success, Adv. Energy Mater., 2021, 11, 2002653.
24. E. Mazzolini, R. A. Pacalaj, Y. Fu, B. R. Patil, R. Patidar, X. Lu, T. M. Watson, J. R. Durrant, Z. Li and N. Gasparini, Pathways to Upscaling Highly Efficient Organic Solar Cells Using Green Solvents: A Study on Device Photophysics in the Transition from Lab-to-Fab, Adv. Sci., 2024, 11, 2402637.
25. B. Ma, Y. Yan, M. Wu, S. Li, M. Ru, Z. Xu and W. Zhao, Sustainable Solution Processing Toward High-Efficiency Organic Solar Cells: A Comprehensive Review of Materials, Strategies, and Applications, Adv. Funct. Mater., 2025, 35, 2413814.
26. 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, Organic solar cells with 20.82% efficiency and high tolerance of active layer thickness through crystallization sequence manipulation, Nat. Mater., 2025, 24, 444–453.
27. M. Azzouzi, J. Yan, T. Kirchartz, K. Liu, J. Wang, H. Wu and J. Nelson, Nonradiative Energy Losses in Bulk-Heterojunction Organic Photovoltaics, Phys. Rev. X, 2018, 8, 031055.
28. X.-K. Chen, D. Qian, Y. Wang, W. Tress, H. Yao, K. Vandewal, O. IInganäs, F. Gao, V. Coropceanu and J.-L. Brédas, A Unified Description of Non-Radiative voltage Losses in Organic Solar Cells, Nat. Energy, 2021, 6, 799–806.
29. Y. Tamai, R. Shirouchi, T. Saito, K. Kohzuki and S.-I. Natsuda, Role of the Energy Offset in the Charge Photogeneration and Voltage Loss of Nonfullerene Acceptor-Based Organic Solar Cells, J. Mater. Chem. A, 2023, 11, 17581–17593.
30. A. Landi, D. Padula and A. Peluso, Fast Nonradiative Decay Paths in Organic Solar Cells: Implications for Designing More Efficient Photovoltaic Systems, ACS Appl. Energy Mater., 2024, 7, 707–714.
31. T. Pizza, A. Landi, F. Ambrosio, A. Capobianco and A. Peluso, Elementary Processes in Ternary Solar Cells, Mater. Chem. Front., 2024, 8, 4069–4076.
32. Y. Cui, P. Zhu, X. Liao and Y. Chen, Recent Advances of Computational Chemistry in Organic Solar Cell Research, J. Mater. Chem. C, 2020, 8, 15920–15939.
33. M. Azzouzi, N. P. Gallop, F. Eisner, J. Yan, X. Zheng, H. Cha, Q. He, Z. Fei, M. Heeney, A. A. Bakulin and J. Nelson, Reconciling Models of Interfacial State Kinetics and Device Performance in Organic Solar Cells: Impact of the Energy Offsets on the Power Conversion Efficiency, Energy Environ. Sci., 2022, 15, 1256–1270.
34. S. M. Menke, N. A. Ran, G. C. Bazan and R. H. Friend, Understanding Energy Loss in Organic Solar Cells: Toward a New Efficiency Regime, Joule, 2018, 2, 23–25.
35. L. Zhu, Z. Tu, Y. Yi and Z. Wei, Achieving Small Exciton Binding Energies in Small Molecule Acceptors for Organic Solar Cells: Effect of Molecular Packing, J. Phys. Chem. Lett., 2019, 10, 4888–4894.
36. X. Zhu, K. Wang, J. He, L. Zhang, H. Yu, D. He and B. Hu, Exploring Deep and Shallow Trap States in a Non-Fullerene Acceptor ITIC-Based Organic Bulk Heterojunction Photovoltaic System, J. Phys. Chem. C, 2019, 123, 20691–20697.
37. D. Qian, Z. Zheng, H. Yao, W. Tress, T. R. Hopper, S. Chen, S. Li, J. Liu, S. Chen, J. Zhang, X.-K. Liu, B. Gao, L. Ouyang, Y. Jin, G. Pozina, I. A. Buyanova, W. M. Chen, O. Inganäs, V. Coropceanu, J.-L. Brédas, H. Yan, J. Hou, F. Zhang, A. A. Bakulin and F. Gao, Design Rules for Minimizing Voltage Losses in High-Efficiency Organic Solar cells, Nat. Mater., 2018, 17, 703–709.
38. N. Zarrabi, O. J. Sandberg, P. Meredith and A. Armin, Subgap Absorption in Organic Semiconductors, J. Phys. Chem. Lett., 2023, 14, 3174–3185.
39. J. Benduhn, K. Tvingstedt, F. Piersimoni, S. Ullbrich, Y. Fan, M. Tropiano, K. A. McGarry, O. Zeika, M. K. Riede, C. J. Douglas, S. Barlow, S. R. Marder, D. Neher, D. Spoltore and K. Vandewal, Intrinsic Non-Rdiative Voltage Losses in Fullerene-Based Organic Solar Cells, Nat. Energy, 2017, 2, 17053.
40. M. Azzouzi, T. Kirchartz and J. Nelson, Factors Controlling Open-Circuit Voltage Losses in Organic Solar Cells, Trends Chem., 2019, 1, 49–62.
41. Y. Miyake, K. Kranthiraja, F. Ishiwari and A. Saeki, Improved Predictions of Organic Photovoltaic Performance through Machine Learning Models Empowered by Artificially Generated Failure Data, Chem. Mater., 2022, 34, 6912–6920.
42. M. Hußner, R. A. Pacalaj, G. O. Müller-Dieckert, C. Liu, Z. Zhou, N. Majeed, S. Greedy, I. Ramirez, N. Li, S. M. Hosseini, C. Uhrich, C. J. Brabec, J. R. Durrant, C. Deibel and R. C. I. MacKenzie, Machine Learning for Ultra High Throughput Screening of Organic Solar Cells: Solving the Needle in the Haystack Problem, Adv. Energy Mater., 2024, 14, 2303000.
43. A. Eibeck, D. Nurkowski, A. Menon, J. B. An Jinkui Wu, L. Zhou, S. Mosbach, J. Akroyd and M. Kraft, Predicting Power Conversion Efficiency of Organic Photovoltaics: Models and Data Analysis, ACS Omega, 2014, 6, 23764–23775.
44. H. Wang, J. Feng, Z. Dong, L. Jin, M. Li, J. Yuan and Y. Li, Efficient Screening Framework for Organic Solar Cells with Deep Learning and Ensemble Learning, npj Comput. Mater., 2023, 9, 200.
45. J. Yan, X. Rodríguez-Martínez, D. Pearce, H. Douglas, D. Bili, M. Azzouzi, F. Eisner, A. Virbule, E. Rezasoltani, V. Belova, B. Dörling, S. Few, A. A. Szumska, X. Hou, G. Zhang, H.-L. Yip, M. Campoy-Quiles and J. Nelson, Identifying Structure-Absorption Relationships and Predicting Absorption Strength of Non-Fullerene Acceptors for Organic Photovoltaics, Energy Environ. Sci., 2022, 15, 2958–2973.
46. O. H. Omar, M. del Cueto, T. Nematiaram and A. Troisi, High-throughput virtual screening for organic electronics: a comparative study of alternative strategies, J. Mater. Chem. C, 2021, 9, 13557–13583.
47. M. Seifrid, S. Lo, D. G. Choi, G. Tom, M. L. Le, K. Li, R. Sankar, H.-T. Vuong, H. Wakidi, A. Yi, Z. Zhu, N. Schopp, A. Peng, B. R. Luginbuhl, T.-Q. Nguyen and A. Aspuru-Guzik, Beyond Molecular Structure: Critically Assessing Machine Learning for Designing Organic Photovoltaic Materials and Devices, J. Mater. Chem. A, 2024, 12, 14540–14558.
48. W. Tress, N. Marinova, O. Inganäs, M. K. Nazeeruddin, S. M. Zakeeruddin and M. Graetzel, Predicting the Open-Circuit Voltage of CH₃NH₃PbI₃ Perovskite Solar Cells Using Electroluminescence and Photovoltaic Quantum Efficiency Spectra: The Role of Radiative and Non-Radiative Recombination, Adv. Energy Mater., 2015, 5, 1400812.
49. R. Borrelli and A. Peluso, The Temperature Dependence of Radiationless Transition Rates from Ab Initio Computations, Phys. Chem. Chem. Phys., 2011, 13, 4420–4426.
50. R. Borrelli, A. Capobianco and A. Peluso, Hole Hopping Rates in Single Strand Oligonucleotides, Chem. Phys., 2014, 440, 25–30.
51. G. Chen, C. Si, Z. Tang, K. Guo, T. Wang, J. Zhang and B. Wei, Temperature-Dependent Device Performance of Organic Photovoltaic Cells Based on a Squaraine Dye, Synth. Met., 2016, 222, 293–298.
52. Q. Burlingame, G. Zanotti, L. Ciammaruchi, E. A. Katz and S. R. Forrest, Outdoor Operation of Small-molecule Organic Photovoltaics, Org. Electron., 2017, 41, 274–279.
53. J. Yao, T. Kirchartz, M. S. Vezie, M. A. Faist, W. Gong, Z. He, H. Wu, J. Troughton, T. Watson, D. Bryant and J. Nelson, Quantifying Losses in Open-Circuit Voltage in Solution-Processable Solar Cells, Phys. Rev. Appl., 2015, 4, 014020.
54. W. Shockley and H. J. Queisser, Detailed Balance Limit of Efficiency of p–n Junction Solar Cells, J. Appl. Phys., 1961, 32, 510–519.
55. J. F. Guillemoles, T. Kirchartz, D. Cahen and U. Rau, Guide for the Perplexed to the Shockley-Queisser Model for Solar Cells, Nat. Photonics, 2019, 13, 501–505.
56. O. J. Sandberg, A. Sundqvist, M. Nyman and R. Österbacka, Relating Charge Transport, Contact Properties, and Recombination to Open-Circuit Voltage in Sandwich-Type Thin-Film Solar Cells, Phys. Rev. Appl., 2016, 5, 044005.
57. T. Kirchartz, F. Deledalle, P. S. Tuladhar, J. R. Durrant and J. Nelson, On the Differences between Dark and Light Ideality Factor in Polymer: Fullerene Solar Cells, J. Phys. Chem. Lett., 2013, 4, 2371–2376.
58. U. Rau, Reciprocity Relation between Photovoltaic Quantum Efficiency and Electroluminescent Emission of Solar Cells, Phys. Rev. B: Condens. Matter Mater. Phys., 2007, 76, 085303.
59. F. A. Lindholm, J. G. Fossum and E. L. Burgess, Application of the Superposition Principle to Solar-Cell Analysis, IEEE Trans. Electron Devices, 1979, 26, 165–171.
60. S. J. Robinson, A. G. Aberle and M. A. Green, Departures from the Principle of Superposition in Silicon Solar Cells, J. Appl. Phys., 1994, 76, 7920–7930.
61. B. Blank, T. Kirchartz, S. Lany and U. Rau, Selection Metric for Photovoltaic Materials Screening Based on Detailed-Balance Analysis, Phys. Rev. Appl., 2017, 8, 024032.
62. National Renewable Energy Laboratory (NREL), AM1.5, https://www.nrel.gov/grid/solar-resource/spectra-am1.5.html, (accessed 2026-04-01).
63. C. A. Gueymard, D. Myers and K. Emery, Proposed Reference Irradiance Spectra for Solar Energy Systems Testing, Sol. Energy, 2002, 73, 443–467.
64. M. Presselt, F. Herrmann, S. Shokhovetz, H. Hoppe, E. Runge and G. Gobbsh, Sub-Bandgap Absorption in Polymer-Fullerene Solar Cells Studied by Temperature-Dependent External Quantum Efficiency and Absorption Spectroscopy, Chem. Phys. Lett., 2012, 542, 70–73.
65. Y. Tamai, Charge Generation in Organic Solar Cells: Journey toward 20% Power Conversion Efficiency, Aggregate, 2022, 3, e280.
66. J. J. Sakurai and J. Napolitano, Modern Quantum Mechanics, Pearson, 2nd edn, 2017.
67. C. X. Zhao, S. Xiao and G. Xu, Density of Organic Thin Films in Organic Photovoltaics, J. Appl. Phys., 2015, 118, 044510.
68. A. Velardo, R. Borrelli, A. Peluso and A. Capobianco, First-Principle Calculations of the Band Shapes of Singlet-Triplet Transitions, J. Phys. Chem. C, 2016, 120, 24605–24614.
69. A. Velardo, R. Borrelli, A. Capobianco, A. Landi and A. Peluso, Disentangling Electronic and Vibrational Effects in the Prediction of Band Shapes for Singlet-Triplet Transitions, J. Phys. Chem. C, 2019, 123, 14173–14179.
70. A. Landi, A. Landi, A. Leo and A. Peluso, The Rates of Non-Adiabatic Processes in Large Molecular Systems: Toward an Effective Full-Dimensional Quantum Mechanical Approach, J. Chem. Phys., 2024, 160, 174114.
71. M. Lax, The Franck–Condon Principle and Its Application to Crystals, J. Chem. Phys., 1952, 20, 1752–1760.
72. R. Kubo and Y. Toyozawa, Application of the Method of Generating Function to Radiative and Non-Radiative Transitions of a Trapped Electron in a Crystal, Prog. Theor. Phys., 1955, 13, 160–182.
73. R. Borrelli, A. Capobianco and A. Peluso, Generating Function Approach to the Calculation of Spectral Band Shapes of Free-Base Chlorin Including Duschinsky and Herzberg-Teller Effects, J. Phys. Chem. A, 2012, 116, 9934–9940.
74. R. Borrelli and A. Peluso, Elementary Electron Transfer Reactions: From Basic Concepts to Recent Computational Advances, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2013, 3, 542–559.
75. T. Pizza, A. Capobianco and A. Troisi, The Importance of High-Frequency Modes in the Prediction of RISC Rates for TADF Molecules, J. Phys. Chem. Lett., 2025, 16, 3056–3062.
76. U. Rau, U. W. Paetzold and T. Kirchartz, Thermodynamics of Light Management in Photovoltaic Devices, Phys. Rev. B: Condens. Matter Mater. Phys., 2014, 90, 035211.
77. W. Wan Roosbroeck and W. Shockley, Photon-Radiative Recombination of Electrons and Holes in Germanium, Phys. Rev., 1954, 94, 1558–1560.
78. Y. Song, B. Li, S. Liu, M. Qin, Y. Gao, K. Zhang, L. Lin, C.-K. Wang and J. Fan, Theoretical Studies on the Excited-State Properties of Thermally Activated Delayed Fluorescence Molecules with Aggregation Induced Emission, J. Mater. Chem. C, 2022, 10, 9377–9390.
79. S. R. Yost, L.-P. Wang and T. Van Voorhis, Molecular Insight Into the Energy Levels at the Organic Donor/Acceptor Interface: A Quantum Mechanics/Molecular Mechanics Study, J. Phys. Chem. C, 2011, 115, 14431–14436.
80. G. Long, A. Li, R. Shi, Y.-C. Zhou, X. Yang, Y. Zuo, W.-R. Wu, U.-S. Jeng, Y. Wang, X. Wan, P. Shen, H.-L. Zhang, T. Yan and Y. Chen, The Evidence for Fullerene Aggregation in High-Performance Small-Molecule Solar Cells by Molecular Dynamics Simulation, Adv. Electron. Mater., 2015, 1, 1500217.
81. A. Troisi and G. Orlandi, The Hole Transfer in DNA: Calculation of Electron Coupling between Close Bases, Chem. Phys. Lett., 2001, 344, 509–518.
82. A. Landi and A. Troisi, Rapid Evaluation of Dynamic Electronic Disorder in Molecular Semiconductors, J. Phys. Chem. C, 2018, 122, 18336–18345.
83. T. Nematiaram, D. Padula, A. Landi and A. Troisi, On the Largest Possible Mobility of Molecular Semiconductors and How to Achieve It, Adv. Funct. Mater., 2020, 30, 2001906.
84. 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. Bredas, H.-L. Yip and Y. Cao, Delocalization of Exciton and Electron Wavefunction in non-Fullerene Acceptor Molecules Enables Efficient Organic Solar Cells, Nat. Commun., 2020, 11, 3943.
85. D. Padula, A. Landi and G. Prampolini, Assessing alkyl side chain effects on electron transport properties of Y6-derived non-fullerene acceptors, Energy Adv., 2023, 2, 1215–1224.
86. Z.-W. Zhao, A. Troisi, Y. Geng, X.-F. Wang and Z.-M. Su, Exploring Charge Dissociation in a Statistical Sample of Active-Layer Models of an Organic Solar Cell, J. Phys. Chem. C, 2020, 124, 18840–18846.
87. F. G. Guijarro, R. Caballero, P. de la Cruz, R. Singhal, F. Langa and G. D. Sharma, Ternary All-Small-Molecule Solar Cells with Two Small-Molecule Donors and Y6 Nonfullerene Acceptor with a Power Conversion Efficiency over Above 14% Processed from a Nonhalogenated Solvent, Solar RRL, 2020, 4, 2000460.
88. A. K. K. Kyaw, D. H. Wang, D. Wynands, J. Zhang, T.-Q. Nguyen, G. C. Bazan and A. J. Heeger, Improved Light Harvesting and Improved Efficiency by Insertion of an Optical Spacer (ZnO) in Solution-Processed Small-Molecule Solar Cells, Nano Lett., 2013, 13, 3796–3801.
89. J. Ge, L. Xie, R. Peng, B. Fanady, J. Huang, W. Song, T. Yan, W. Zhang and Z. Ge, 13.34% Efficiency Non-Fullerene All-Small-Molecule Organic Solar Cells Enabled by Modulating the Crystallinity of Donors via a Fluorination Strategy, Angew. Chem., Int. Ed., 2020, 59, 2808–2815.
90. B. Kan, M. Li, Q. Zhang, F. Liu, X. Wan, Y. Wang, W. Ni, G. Long, X. Yang, H. Feng, Y. Zuo, M. Zhang, F. Huang, Y. Cao, T. P. Russell and Y. Chen, A Series of Simple Oligomer-like Small Molecules Based on Oligothiophenes for Solution-Processed Solar Cells with High Efficiency, J. Am. Chem. Soc., 2015, 137, 3886–3893.
91. J. Yuan and Y. Zou, The History and Development of Y6, Org. Electron., 2022, 102, 106436.
92. A. Melianas, V. Pranculis, D. Spoltore, J. Benduhn, O. Inganäs, V. Gulbinas, K. Vandewal and M. Kemerink, Charge Transport in Pure and Mixed Phases in Organic Solar Cells, Adv. Energy Mater., 2017, 7, 1700888.
93. S. Kazaoui, N. Minami, Y. Tanabe, H. J. Byrne, A. Eilmes and P. Petelenz, Comprehensive Analysis of Intermolecular Charge-Transfer Excited States in C₆₀ and C₇₀ Films, Phys. Rev. B: Condens. Matter Mater. Phys., 1998, 58, 7689–7700.
94. S. J. Finkelmeyer and M. Presselt, Tuning Optical Properties of Organic Thin Films through Intermolecular Interactions – Fundamentals, Advances and Strategies, Chem. – Eur. J., 2025, 31, e202403500.
95. S. Alam, V. Nádaždy, T. Váry, C. Friebe, R. Meitzner, J. Ahner, A. Anand, S. Karuthedath, C. S. P. De Castro, C. Göhler, S. Dietz, J. Cann, C. Kästner, A. Konkin, W. Beenken, A. M. Anton, C. Ulbricht, A. Sperlich, M. D. Hager, U. Ritter, F. Kremer, O. Brüggemann, U. S. Schubert, D. A. M. Egbe, G. C. Welch, V. Dyakonov, C. Deibel, F. Laquai and H. Hoppe, Uphill and Downhill Charge Generation from Charge Transfer to Charge Separated States in Organic Solar Cells, J. Mater. Chem. C, 2021, 9, 14463–14489.
96. Y. Dong, V. C. Nikolis, F. Talnack, Y.-C. Chin, J. Benduhn, G. Londi, J. Kublitski, X. Zheng, S. C. B. Mannsfeld, D. Spoltore, L. Muccioli, J. Li, X. Blase, D. Beljonne, J.-S. Kim, A. A. Bakulin, G. D'Avino, J. R. Durrant and K. Vandewal, Orientation Dependent Molecular Electrostatics Drives Efficient Charge Generation in Homojunction Organic Solar Cells, Nat. Commun., 2020, 11, 4617.
97. K. Eguchi and H. Murata, The Ionization Energy of α-Sexithiophene and p-Sexiphenyl in 2D and 3D Thin Films grown on Silicon Oxide Surfaces, Phys. Chem. Chem. Phys., 2024, 26, 8687–8694.
98. C. Jehoulet, A. J. Bard and F. Wudl, Electrochemical Reduction and Oxidation of C60 films, J. Am. Chem. Soc., 1991, 113, 5456–5457.
99. M. Namazian, C. Y. Lin and M. L. Coote, Benchmark Calculations of Absolute Reduction Potential of Ferricinium/Ferrocene Couple in Nonaqueous Solutions, J. Chem. Theory Comput., 2010, 6, 2721–2725.
100. M. I. Pilo, E. Masolo, L. Maidich, P. Manca, G. Sanna, N. Spano and A. Zucca, Voltammetric and Spectroscopic Investigation of Electrogenerated Oligo-Thiophenes: Effect of Substituents on the Energy-Gap Value, Appl. Sci., 2022, 12, 11714.
101. G. Zotti, B. Vercelli, A. Berlin, S. Destri, M. Pasini, V. Hernández and J. L. Navarrete, Electrochemical, Magnetic, and Electrical Properties of α,ω-Capped Sexithiophene Films. Part 3. Conduction in Poly(bis-terthienyl-B)s (B = Ethane, Disulfide, Diacetylene, Acetylene, Ethylene), Chem. Mater., 2008, 20, 6847–6856.
102. J. Zhang, H. S. Tan, X. Guo, A. Facchetti and H. Yan, Material Insights and Challenges for Non-Fullerene Organic Solar Cells Based on Small Molecular Acceptors, Nat. Energy, 2018, 3, 720–731.
103. M. Jeong, S. Chen, S. M. Lee, Z. Wang, Y. Yang, Z.-G. Zhang, C. Zhang, M. Xiao, Y. Li and C. Yang, Feasible D1–A–D2–A Random Copolymers for Simultaneous High-Performance Fullerene and Nonfullerene Solar Cells, Adv. Energy Mater., 2018, 8, 1702166.
104. Q. Zhang, X. Yuan, Y. Feng, B. W. Larson, G. M. Su, Y. Maung Maung, N. Rujisamphan, Y. Li, J. Yuan and W. Ma, Understanding the Interplay of Transport-Morphology-Performance in PBDB-T-Based Polymer Solar Cells, Solar RRL, 2020, 4, 1900524.
105. P. E. Schwenn, P. L. Burn and B. J. Powell, Calculation of Solid State Molecular Ionisation Energies and Electron Affinities for Organic Semiconductors, Org. Electron., 2011, 12, 394–403.
106. A. Capobianco, J. Wiktor, A. Landi, F. Ambrosio and A. Peluso, Electron Localization and Mobility in Monolayer Fullerene Networks, Nano Lett., 2024, 24, 8335–8342.
107. M. Yanxia, Q. L. Liang, H. W. Wu and Y. C. Cao, High Photogeneration and Low Recombination Rate Leading to High-Performance Non-Fullerene Organic Solar Cells, J. Mater. Chem. A, 2023, 11, 6237–6247.
108. N. Wei, Y. Guo, H. Song, Y. Liu, H. Lu and Z. Bo, Reducing Non-Radiative Energy Losses in Non-Fullerene Organic Solar Cells, ChemSusChem, 2025, 18, e202402169.
109. J. Hofinger, C. Putz, F. Mayr, K. Gugujonovic, D. Wielend and M. C. Scharber, Understanding the Low Voltage Losses in High-Performance Non-Fullerene Acceptor-Based Organic Solar Cells, Mater. Adv., 2021, 2, 4291–4302.
110. R. Englman and J. Jortner, The Energy Gap Law for Radiationless Transitions in Large Molecules, Mol. Phys., 1970, 18, 145–164.
111. I. R. Gould, D. Noukakis, L. Gomez-Jahn, R. H. Young, J. L. Goodman and S. Farid, Radiative and Nonradiative Electron Transfer in Contact Radical-Ion Pairs, Chem. Phys., 1993, 176, 439–456.
112. U. Hörmann, C. Lorch, A. Hinderhofer, A. Gerlach, M. Gruber, J. Kraus, B. Sykora, S. Grob, T. Linderl and A. Wilke, et al., V_OC from a Morphology Point of View: The Influence of Molecular Orientation on the Open Circuit Voltage of Organic Planar Heterojunction Solar Cells, J. Phys. Chem. C, 2014, 118, 26462–26470.
113. J. Widmer, K. Leo and M. Riede, Temperature Dependent Behavior of Flat and Bulk Heterojunction Organic Solar Cells, MRS Online Proc. Libr., 2013, 1493, 141–145.
114. M. Lenes, L. J. A. Koster, V. D. Mihailetchi and P. W. M. Blom, Thickness Dependence of the Efficiency of Polymer: Fullerene Bulk Heterojunction Solar Cells, Appl. Phys. Lett., 2006, 88, 243502.
115. X. Liu, Y. Li, P. Huang, Y. Zhou, Z.-Q. Jiang, B. Song, Y. Li, L.-S. Liao and Y. Zheng, Highly Efficient and Thickness-Tolerable Bulk Heterojunction Polymer Solar Cells Based on P3HT donor and a Low-Bandgap Non-Fullerene Acceptor, J. Power Sources, 2017, 364, 426–431.
116. C. Sun, F. Pan, H. Bin, J. Zhang, L. Xue, B. Qiu, Z. Wei, Z.-G. Zhang and Y. Li, A Low Cost and High Performance Polymer Donor Material for Polymer Solar Cells, Nat. Commun., 2018, 9, 743.
117. X. Wu, X. Xie and A. Troisi, Calibration of Several First Excited State Properties for Organic Molecules through Systematic Comparison of TDDFT with Experimental Spectra, J. Mater. Chem. C, 2024, 12, 18886–18892.
118. M.-C. Kim, E. Sim and K. Burke, Communication: Avoiding Unbound Anions in Density Functional Calculations, J. Chem. Phys., 2011, 134, 171103.
119. F. Jensen, Describing Anions by Density Functional Theory: Fractional Electron Affinity, J. Chem. Theory Comput., 2010, 6, 2726–2735.
120. J. Simons, Molecular Anions Perspective, J. Phys. Chem. A, 2023, 127, 3940–3957.
121. J. Tomasi, B. Mennucci and R. Cammi, Quantum Mechanical Continuum Solvation Models, Chem. Rev., 2005, 105, 2999–3093.
122. W. Barford and M. Marcus, Theory of Optical Transitions in Conjugated Polymers. I. Ideal Systems, J. Chem. Phys., 2014, 141, 164101.
123. M. Marcus, O. R. Tozer and W. Barford, Theory of Optical Transitions in Conjugated Polymers. II. Real Systems, J. Chem. Phys., 2014, 141, 164102.
124. W. Barford and O. R. Tozer, Theory of Exciton Transfer and Diffusion in Conjugated Polymers, J. Chem. Phys., 2014, 141, 164103.
125. M. A. Green, Lambertian Light Trapping in Textured Solar Cells and Light-Emitting Diodes: Analytical Solutions, Prog. Photovoltics, 2002, 10, 235–241.
126. S. Wilken, T. Upreti, A. Melianas, S. Dahlström, G. Persson, E. Olsson, R. Österbacka and M. Kemerink, Experimentally Calibrated Kinetic Monte Carlo Model Reproduces Organic Solar Cell Current-Voltage Curve, Solar RRL, 2020, 4, 2000029.
127. W. Kaiser, J. Popp, M. Rinderle, T. Albes and A. Gagliardi, Generalized Kinetic Monte Carlo Framework for Organic Electronics, Algorithms, 2018, 11, 37.
128. D. Balzer and I. Kassal, Mechanism of Delocalization-Enhanced Exciton Transport in Disordered Organic Semiconductors, J. Phys. Chem. Lett., 2023, 14, 2155–2162.
129. M. J. Abraham, T. Murtola, R. Schulz, S. Páll, J. C. Smith, B. Hess and E. Lindahl, GROMACS: high performance molecular simulations through multi-level parallelism from laptops to supercomputers, SoftwareX, 2015, 1–2, 19–25.
130. W. L. Jorgensen, D. S. Maxwell and J. Tirado-Rives, Development and Testing of the OPLS All-Atom Force Field on Conformational Energetics and Properties of Organic Liquids, J. Am. Chem. Soc., 1996, 118, 11225–11236.
131. T. Darden, D. York and L. Pedersen, Particle mesh Ewald: an N·log(N) method for Ewald sums in large systems, J. Chem. Phys., 1993, 98, 10089–10092.
132. B. Hess, H. Bekker, H. J. C. Berendsen and J. G. E. M. Fraaije, LINCS: a linear constraint solver for molecular simulations, J. Comput. Chem., 1997, 18, 1463–1472.
133. H. J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren, A. DiNola and J. R. Haak, Molecular dynamics with coupling to an external bath, J. Chem. Phys., 1984, 81, 3684–3690.
134. M. Parrinello and A. Rahman, Polymorphic transitions in single crystals: a new molecular dynamics method, J. Appl. Phys., 1981, 52, 7182–7190.
135. A. D. Becke, Density-Functional Thermochemistry. III. The Role of Exact Exchange, J. Chem. Phys., 1993, 98, 5648–5652.
136. P. J. Stephens, F. J. Devlin, C. F. Chabalowski and M. J. Frisch, Ab Initio Calculation of Vibrational Absorption and Circular Dichroism Spectra Using Density Functional Force Fields, J. Phys. Chem., 1994, 98, 11623–11627.
137. G. D'Avino, L. Muccioli, C. Zannoni, D. Beljonne and Z. G. Soos, Electronic Polarization in Organic Crystals: A Comparative Study of Induced Dipoles and Intramolecular Charge Redistribution Schemes, J. Chem. Theory Comput., 2014, 10, 4959–4971.
138. A. Landi, R. Borrelli, A. Capobianco, A. Velardo and A. Peluso, Hole Hopping Rates in Organic Semiconductors: A Second-Order Cumulant Approach, J. Chem. Theory Comput., 2018, 14, 1594–1601.
139. A. Dreuw and M. Head-Gordon, Single-Reference ab Initio Methods for the Calculation of Excited States of Large Molecules, Chem. Rev., 2005, 105, 4009–4037.
140. M. J. G. Peach, P. Benfield, T. Helgaker and D. J. Tozer, Excitation Energies in Density Functional Theory: An Evaluation and a Diagnostic Test, J. Chem. Phys., 2008, 128, 44118.
141. R. J. Magyar and S. Tretiak, Dependence of Spurious Charge-Transfer Excited States on Orbital Exchange in TDDFT: Large Molecules and Clusters, J. Chem. Theory Comput., 2007, 3, 976–987.
142. D. Mester and M. Kállay, Charge-Transfer Excitations within Density Functional Theory: How Accurate Are the Most Recommended Approaches?, J. Chem. Theory Comput., 2022, 18, 1646–1662.
143. D. Mester and M. Kállay, Accurate Spectral Properties within Double-Hybrid Density Functional Theory: A Spin-Scaled Range-Separated Second-Order Algebraic-Diagrammatic Construction-Based Approach, J. Chem. Theory Comput., 2022, 18, 865–882.
144. T. Yanai, D. P. Tew and N. C. Handy, A New Hybrid Exchange-Correlation Functional Using the Coulomb-Attenuating Method (CAM-B3LYP), Chem. Phys. Lett., 2004, 393, 51–57.
145. P. Dev, S. Agrawal and N. J. English, Determining the Appropriate Exchange-Correlation Functional for Time-Dependent Density Functional Theory Studies of Charge-Transfer Excitations in Organic Dyes, J. Chem. Phys., 2012, 136, 224301.
146. A. Mandal, C. E. Mohn, C. H. Görbitz and A. Roy, Cocrystal Engineering of Organic Semiconductors for Photovoltaic Applications: Modeling Excited-State Properties of a Charge Transfer Cocrystal of a Dicarbazole Donor and a Fluoranil Acceptor, J. Phys. Chem. C, 2026, 130, 882–897.
147. S. Giannini, W. Peng, L. Cupellini, D. Padula, A. Carof and J. Blumberger, Exciton Transport in Molecular Organic Semiconductors Boosted by Transient Quantum Delocalization, Nat. Commun., 2022, 13, 2755.
148. M. B. Price, P. A. Hume, A. Ilina, I. Wagner, R. R. Tamming, K. E. Thorn, W. Jiao, A. Goldingay, P. J. Conaghan, G. Lakhwani, N. J. L. K. Davis, Y. Wang, P. Xue, H. Lu, K. Chen, X. Zhan and J. M. Hodgkiss, Free Charge Photogeneration in a Single Component High Photovoltaic Efficiency Organic Semiconductor, Nat. Commun., 2022, 13, 2827.
149. 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.
150. S. Grimme, J. Antony, S. Ehrlich and H. Krieg, A Consistent and Accurate Ab Initio Parametrization of Density Functional Dispersion Correction (DFT-D) for the 94 Elements H-Pu, J. Chem. Phys., 2010, 132, 154104–154119.
151. J. Aragó and A. Troisi, Excitonic Couplings between Molecular Crystal Pairs by a Multistate Approximation, J. Chem. Phys., 2015, 142, 164107.
152. K. A. Kistler, F. C. Spano and S. Matsika, A Benchmark of Excitonic Couplings Derived from Atomic Transition Charges, J. Phys. Chem. B, 2013, 117, 2032–2044.
153. P. H. Schönemann, A Generalized Solution of the Orthogonal Procrustes Problem, Psychometrika, 1966, 31, 1–10.
154. F. Duschinsky, Zur Deutung der Elektronenspektren Mehratomiger Moleküle. I. Uber das Franck–Condon Prinzip, Acta Physicochim. URSS, 1937, 7, 551–566.
155. R. Borrelli, A. Capobianco and A. Peluso, MolFC: A program for Franck–Condon integrals calculation, Package available online at https://github.com/rborrelli/molfc.
156. A. Capobianco, R. Borrelli, C. Noce and A. Peluso, Franck–Condon Factors in Curvilinear Coordinates: The Photoelectron Spectrum of Ammonia, Theor. Chem. Acc., 2012, 131, 1181.
157. A. Peluso, R. Borrelli and A. Capobianco, Photoelectron Spectrum of Ammonia, a Test Case for the Calculation of Franck–Condon Factors in Molecules Undergoing Large Geometrical Displacements upon Photoionization, J. Phys. Chem. A, 2009, 113, 14831–14837.

---