Effect of thermal disorder on the electronic structure and the charge mobility of acenes

Abstract

The effect of thermal disorder on the electronic structure and transport properties of two prototypical organic semiconductors, naphthalene and pentacene, has been studied by combining molecular dynamics simulations with ab initio electronic structure calculations. Our approach explicitly account for the influence of thermal disorder on the site energy and on the intermolecular interactions. It is shown that thermal disorder shifts the polaronic energy levels with respect to the perfectly ordered crystal and stabilize charge localization, even in the case of pentacene, with a significant contribution provided by local rearrangement of the first neighbours. Evaluation of charge mobilities, carried out by kinetic Monte Carlo simulations, shows that site energy fluctuations play the most significant role in slowing down charge diffusion, leading to computed mobilities which are in very good agreement with the experimental ones.

---

Introduction

Organic semiconductors (OS) based on π-conjugated systems are low-cost alternatives to inorganic materials, presenting several advantages over them, such as easy manufacturing, the possibility of finely tuning their physico-chemical properties for the application of interest, bio-compatibility and flexibility.^1–9 However, the main drawback of OS based devices is their lower efficiencies with respect to their inorganic counterparts.^10,11

Two key features characterize charge transport in organic semiconductors: (i) low reorganization energies due to small dielectric constants and (ii) large amplitude thermal motions, which could significantly modulate charge mobility, by causing time fluctuations of both molecular energy levels and the electronic couplings between them.^12,13

These peculiarities have led to uncertainties regarding the actual mechanism of charge transport in OSs. In fact, even though several theoretical approaches have been able to predict results in good agreement with the experimental charge mobility,^14–19 neither the band-like mechanism nor the hopping mechanism appear to be fully satisfying.^{12,16,20–25} Indeed, while the former has been deemed unsuitable due to the short mean free path observed in organic field-effect transistors,¹³ the latter often overestimates charge mobilities and, when combined with rate constants derived from classical Marcus theory, predicts a temperature dependence that doesn't properly align with the observed trends.²⁶ However, the correct temperature dependence can be recovered by properly incorporating quantum mechanical effects into the calculations of the hopping rates, leading to mobilities which are no longer thermally activated,^15,27,28 thus indicating that the temperature dependence of mobility alone cannot be used to discriminate between band-like and hopping transport mechanisms. Moreover, mobility overestimation could be accommodated by properly accounting for the effects of structural disorder, due to defects and to large amplitude thermal motion. Indeed, the fluctuations of the transfer integral, induced by thermal motion of the molecules around their equilibrium positions, have been considered the main factor limiting the charge mobility of OSs,^1,12 while, in this context, the effect of thermal disorder on site energies have received relatively less attention.²⁹ Moreover, the polaron binding energy, i.e. the stabilization of the localized charge upon reorganization of the molecule bearing it and of the surrounding medium, is also influenced by thermal disorder. In fact, not only a distribution of the polaronic energy levels is expected but also the average binding energy can be sensitively affected by disorder, which usually provides further stabilization of the localized state with respect to the ordered material.³⁰

Herein, with the aim of better quantifying the effect of thermal disorder and of the environment reorganization, for the first time we combine molecular dynamics simulations based on quantum-mechanically-derived force fields^31,32 of neutral and charged molecules with advanced electronic-structure calculations to calculate (i) polaron binding energies for localized holes at room-temperature (ii) and charge mobility of naphthalene and pentacene (see Fig. 1), two prototypical systems whose electronic properties have been widely studied in the literature.^{16,18,33–39} Among the stable condensed acene molecules, they are those with the largest (naphthalene) and the smallest (pentacene) reorganization energies,⁴⁰ and, therefore, represent the limit cases for both charge localization and mobility. Our approach allows for a definition of charge localization within a supercell scheme, employing a Koopmans' compliant density functional, which, being essentially devoid of self-interaction error, ensures a correct description of the extent and energetics of charge localization in the materials under study, taking simultaneously into account the local rearrangement of the molecule hosting the extra charge, the reorganization of the surrounding acene molecules and the room-temperature disorder.

Fig. 1 Structural representation of herringbone structure of (a) naphthalene and (b) pentacene. Unique hole hopping paths (differing by symmetry) are reported as blue arrows and identified by letter A, B and C.

Although these effects have been addressed in some studies on charge mobility, they are often considered in a more averaged or simplified manner, without capturing the full precision and detail emphasized in our approach. Even more advanced theoretical frameworks, such as the transient localization theory,^1,13,41 still approximately treat the energy distribution via the top of the valence band and hence do not take into account the actual polaron binding energy. The approach adopted here also allows for a deeper understanding of the effects which concur to determine charge localization in apolar environments, such as those of typical OSs.

Indeed, in this work we will show that thermal disorder shifts the polaronic energy levels with respect to the perfectly ordered crystal, stabilizing charge localization, and that the reorganization of the environment around the charged molecule, evidenced by a contraction of the first neighbours towards the charge, significantly contribute to polaron stabilization. Those effects cannot be captured by implicit solvation models and molecular dynamics simulations based on classical force fields of neutral molecules. In fact, most of the previous studies predict a vanishingly small contribution of the environment to the global reorganization energy,^42,43 whereas more recent theoretical works suggest higher values of solvent reorganization energies in apolar media.^44–46

Finally, we consider both energy and coupling fluctuations obtained by electronic computations to feed a Kinetic Monte Carlo approach for the evaluation of charge mobility,^47–49 where charge transfer rates are evaluated on-the-fly at full quantum mechanical level (using Fermi's Golden Rule). The above approach shows that mobility decay with respect to perfect crystal is primarily due to site energy fluctuations, while electronic coupling fluctuations appear to have a minor effect.

Computational details

The structures of the molecules studied in this work have been retrieved from the Cambridge Structural Database⁵⁰ (CSD); 600187 and 1230799 are the CSD codes for naphthalene and pentacene, respectively.

Classical molecular dynamics

Since in a related paper⁴⁰ we have shown that classical molecular dynamics (MD) with QM-derived force fields lead to results in excellent agreement with much more expensive DFT-based molecular dynamics (DFT-MD) simulations, in this work we have decided to alleviate the computational burden by performing single-point DFT calculations on top of snapshots extracted by MD simulations.

We have first produced accurate intramolecular quantum-mechanically-derived force field (QMD-FF) parameterization for all the molecules under study (either in their neutral and in their cationic form), adopting the JOYCE parameterization protocol,^31,32,51 which exploits a fitting procedure of the Hessian matrix obtained at the same level of theory described in the Section describing the Electronic couplings. In Section S1 of the SI, an extensive description of the parameterization procedure has been reported. The description of the molecules for classical MD simulations has been completed by combining the intramolecular QMD-FFs with Lennard-Jones parameters from appropriate OPLS atom types,⁵² and point charges fitted following the RESP procedure.⁵³ The resulting force fields have been made available to the community in a public GitHub repository (https://github.com/AlessandroUniSa1/acenes_paper2025).

Using these QMD-FFs, we ran classical MD simulations with the 2020.5 software.⁵⁴ We generated a cubic 2 × 2 × 2 box, adopting periodic boundary conditions and taking into account long range electrostatic effects through the PME algorithm.⁵⁵ We used an integration time step of 2 fs, imposing constraints on the bonds involving H atoms through the LINCS algorithm,⁵⁶ adopting a modified Berendsen thermostat⁵⁷ to control temperature and a Parrinello-Rahman barostat to control pressure.⁵⁸ The computational protocol consisted of an initial steepest descent minimisation, followed by a slow NVT heating procedure starting from 25 K, initially carried out keeping aromatic cores frozen, and gradually releasing the restrictions at 300 K. After equilibrating the system at the target temperature in NVT conditions for 0.5 ns, we carried out a production run in NVT conditions, for 2 ns.

In order to compute the quantities in eqn (6), we also need to perform MD simulations for the systems bearing a positive charge. To tackle this task, the geometry of an acene molecule in the center of the 2 × 2 × 2 box used for the neutral system has been replaced with the geometry of the corresponding acene in its cationic form optimized at the B3LYP/6-31G* level. The charged molecule is then described with its own QMD-FF in MD computations. Except from this, the same MD protocol described above for the neutral system has been followed.

We noticed that, while classical MD simulations with larger supercells would imply an affordable computational cost, we limited our analysis to 2 × 2 × 2 supercells because of the high computational overhead associated with hybrid DFT calculations to be performed on top of the MD-produced structural configurations. However, we have performed additional MD simulations for the charged naphthalene supercell, considering a 4 × 4 × 4 supercell and no differences have been noticed with respect to 2 × 2 × 2 supercell, see Section 2.2 in the SI.

Density functional theory calculations on periodic supercells

We perform DFT calculations using a Koopmans' compliant formulation of hybrid functionals: in particular, we use the scheme based on the PBE0 family of functionals.^59,60 The fraction of Fock exchange fulfilling the generalized Koopmans' condition for naphthalene and pentacene have been found to be 37% and 33%, respectively, in ref. 40. We note that this belongs to a class of functionals that have been constructed^61–64 to silence the self-interaction error that affects standard semi-local but also hybrid DFT.^65,66 Therefore, the functional correctly describes the extent and the energetics of charge localization in condensed phase, as previously observed for a plethora of semiconductors and insulators.^61–64 Calculations includes non local-electron correlation via the self-consistent rVV10 scheme.^67,68

Electronic structure calculations on the same supercells employed for classical MD simulations are performed using the freely available CP2K\QUICKSTEP package,⁶⁹ in which core electrons are described by analytical Goedecker-Teter-Hutter pseudopotentials,⁷⁰ while atomic basis sets are employed for valence electrons. In particular, we use the MOLOPT double-zeta polarized basis set⁷¹ and a cutoff of 600 Ry for the plane waves. We take advantage of the auxiliary matrix method, available in CP2K, which enables faster hybrid functional calculations.^72,73 We employ the cFIT3 basis set^72,73 as auxiliary basis to evaluate the two-electron integrals. For both naphthalene and pentacene, the evaluation of 〈ΔE_ox〉_R0 is carried out as follows: we consider 50 structural configurations equally spaced in time from molecular dynamics simulations of the neutral supercell and we optimize the wavefunction for q = 0 and q = + 1; for the latter, we employ the unrestricted Kohn–Sham formalism. Analogously, 〈ΔE_red〉_R+1 is estimated with the same procedure employing structural configurations of the polaron classical MD simulation (i.e. the MD simulation with one charged molecule in the system).

Electronic calculations for electronic couplings

Since the self-consistent charge density-functional tight-binding (DFTB), which we will use to compute supercell phonons, vide infra, is known to severely underestimate electronic couplings because of the minimal basis set used,^16,74 equilibrium geometries, normal modes, and vibrational frequencies of naphthalene and pentacene in their neutral and cationic forms were obtained at the DFT level using the B3LYP functional with the 6-31+G(d,p) basis set, which has been proved to yield sufficiently accurate results.^16,75,76 Dielectric effects have been estimated using the polarizable continuum model (PCM).⁷⁷ According to previous works, we have assumed ε = 4 as the average value for the dielectric constants of the bulk materials.⁷⁸ The Gaussian package has been used for electronic wavefunction computations.⁷⁹

Evaluation of non-local electron-phonon coupling

Following previous works,^16,80,81 we have represented the mode M as a vector of Cartesian displacements Q_M = {x^M_k}, so that S_ij,M can be expressed as:

S_ij,M = ∇J_ij·Q_M. (1)

here, ∇J_ij is the cartesian gradient of the electronic coupling, i.e. the derivative of the electronic coupling with respect to the cartesian displacement of an atom k:where the elements of the vector ∇J_ij are zero if displacements x_k do not belong to the molecules i or j. This approach is particularly suited for the development of a fast protocol, because the first term of the product (∇J_ij) only includes 6N_A differentiations (where N_A is the number of atoms in one molecule) and the second term does not depend on the electronic coupling, so it needs to be evaluated only once.

Supercell phonons

Supercell phonons have been calculated by employing the DFTB method as implemented in the DFTB+ software package.^82,83 The 3-ob-1 Slater-Koster set parameters have been employed.⁸⁴ The atomic positions of the unit cell have been optimized keeping the experimental lattice values fixed, using the conjugate gradient method with a force threshold criterion of 10⁻⁸ Hartree per Bohr radius. We have used periodic boundary conditions employing a 2 × 2 × 2 Monkhorst–Pack k-point sampling scheme,⁸⁵ as suggested in previous works.^37,80,86 The Lennard-Jones dispersion interactions⁸⁷ have been included to take into account the weak van der Waals interaction in the crystalline phase. To investigate the phonon dispersion relations in crystalline phases, we have employed the frozen phonon calculation in which the sampling of the q-space relies on calculation of second derivatives of the energy on sufficiently large supercells. As suggested in previous works,^80,88 a 2 × 2 × 2 supercell has been used. These latter calculations have been carried out employing a 1 × 1 × 1 Monkhorst–Pack k-point sampling scheme.

Charge transfer rates and mobilities

Franck–Condon weighted densities of states for FGR rates' calculations have been computed by using a development version of the MolFC package⁸⁹ available upon request from the authors. The curvilinear coordinate representation of the normal modes has been adopted to prevent that large displacements of an angular coordinate could reflect into large shifts of the equilibrium positions of the involved bond distances. That unphysical effect is unavoidable when using rectilinear coordinates and requires the use of high order anharmonic potentials for its correction.^90–93

As discussed more into details in the SI, the rates used here share the same formal origin as those implemented in MOMAP⁹⁴ by Shuai's group: both are derived from Fermi's Golden Rule and evaluate a vibrational (Franck–Condon) correlation function in the displaced-oscillator (Duschinsky), but they differ from an implementative point of view; thus, where numerical differences occur, they can be traced to the implementation differences.^27,28,94

KMC simulations have been performed using an in-house developed code⁹⁵ similar to the algorithm described in ref. 96, but (i) the computations of the rates of each step have been evaluated by eqn (15) and (ii) only charge transfer processes have been considered.⁹⁷

When included, the influence of disorder on site energies and electronic couplings is accounted for by introducing stochastic variations at the start of each trajectory. In our KMC framework, site energies and electronic couplings are stochastically sampled from Gaussian distributions at the beginning of each trajectory, with parameters derived from hybrid-DFT analyses of MD snapshots (which by their very nature encapsulate dynamic thermal fluctuations in the molecular ensemble). These values remain fixed throughout individual trajectories to model static disorder per realization, but the overall mobility is obtained by averaging over 20 000 such trajectories, each with a different realization of the disorder, providing an ensemble-averaged result that approximates the effects of dynamic disorder via the ergodic principle.^98,99 This method aligns with established approaches in the field, where static ensemble averaging has been demonstrated to yield mobilities consistent with experimental observations in disordered systems, though it may not fully resolve ultra-fast dynamic contributions that could hinder transport in high-mobility regimes.^98,99 To overcome this limitation, future work incorporating explicit delocalization and dynamic effects in time-dependent KMC are ongoing in our group.

We note that our protocol is restricted to keep the charge on one molecule, based on the results of our hybrid-DFT relaxations of representative snapshots, which confirm preserved localization, and IPR analysis which shows that naphthalene's (and the vast majority of pentacene's) snapshots show a molecular localization of the charge. In any case, using a localized basis is a formally valid representation (when the FGR/perturbative assumptions and the relevant timescale separation hold) since an ensemble of fast transitions between localized states would reproduce the macroscopic transport that would be obtained from a delocalized (large-polaron or band) basis. In light of these aspects, we expect that (computationally intensive) explicit two-molecule delocalization or delocalized-KMC could refine rates but are unlikely to alter our qualitative conclusions.

Results and discussion

Hole polarons in disordered acene crystals at room temperature

We first compare the electronic structure of the neutral and positively charged systems to preliminary assess charge localization upon hole injection. To this end, we inspect the respective electronic density of states, as achieved from an average over 50 MD structural configurations equi-spaced in time. Since in a related paper⁴⁰ we have shown that classical molecular dynamics (MD) with QM-derived force fields (QMD-FF) leads to results in excellent agreement with much more expensive DFT-based molecular dynamics simulations, in this work we have decided to alleviate the computational burden by performing single-point DFT calculations on top of snapshots extracted by MD simulations, of which the details are reported in the Computational Details section and in Section S1 in the SI. We here underline that, within the nature of our classical MD framework, single molecule localization is fixed in classical simulations of charged supercells and one could argue that this may affect the subsequent DFT analysis. In this regard, we note that, for naphthalene, hole localization on a single molecule has been demonstrated by ab initio calculations of ref. 40. For pentacene, previous results suggest that localized and delocalized states in the perfect crystal are very close in energy.⁴⁰ For this reason, in this case, to exclude that the occurrence of localization derives from the starting assumption imposed by the force field, we not only reevaluate the total energy of the structural configurations achieved with classical MD but, as an additional test, we perform hybrid-DFT geometry optimizations on top of them. The results, cf. SI, indicate that the degree of localization observed from MD configurations is preserved upon relaxation and is achieved only by a correct treatment of exact exchange in hybrid-DFT. This also suggests that room-temperature disorder brings to the occurrence of energy barriers preventing delocalization of the charge at 0 K, as observed in ref. 40 and that any bias arising from our starting classical MD model should be at least partially mitigated.

For naphthalene, cf. Fig. 2(a), we note the appearance of a second small peak in the DOS of the charged system, centered at ≈0.6 eV above the valence band maximum (VBM) of the material: this clearly denotes a localized state in the gap, i.e. the formation of a hole polaron. In contrast, for the positively charged pentacene, cf. Fig. 2(b), we observe the occurrence of a wide shoulder in the valence band (VB) DOS extending up to ≈0.4 eV above the VBM of the semiconductor. This feature in the DOS presents a small peak closer to the VBM and a broader part at higher energies, suggesting the coexistence of semi-localized and more localized states. We note that, for the ordered material, no charge localization was envisaged in ref. 40, thus suggesting that disorder plays a key role in favouring polaron formation in the close competition with delocalized charge carriers.

Fig. 2 Electronic DOS near the valence band edge for neutral and positively charged (a) naphthalene and (b) pentacene supercells, as achieved from the average of 50 MD structural configurations. Dashed lines point to the value of VBM as linearly extrapolated from the wing of the valence band DOS. Neutral and charged calculations are aligned via the C 2s core level. In each panel, energies are referred to that of the VBM. In panel (b) a darker shade is used to highlight the DOS of the lowest unoccupied molecular orbital.

To delve deeper into this issue, we calculate the inverse participation ratio (IPR) of the wave function of the hole, which provides a quantitative description of how spread out the extra charge is in the supercell:¹⁰⁰

where N is the number of volume elements in the supercell. From eqn (3), if the charge is spread out on a large number of volume elements in the supercell, all |ψ(i)|² terms are small and further squaring them gives even smaller values. In contrast, if there are a few sites on which the amplitude of the wavefunction is sizable, the magnified contributions on |ψ(i)|⁴ leads to larger IPRs. Therefore, within this formulation, an ideally delocalized state should give IPR = 1, while larger values are expected for a localized charge. We then investigate the correlation between hole localization and position of the associated Kohn–Sham (KS) energy level in the band gap, for each structural configuration. To ensure a meaningful comparison, we also calculate the IPR associated with (i) a fully delocalized hole on the VBM and (ii) a hole entirely localized on a single acene molecule. The former is estimated from the hole wave function achieved upon vertical injection onto the ordered material. The latter is determined by performing a calculation with a single acene molecule within the same supercell employed for the condensed-phase system. In this way, we obtain two limit values for IPR. Further, in order to verify whether disorder already localizes to some degree the VB edge of the acene crystal, we also calculate the average IPR of the hole wavefunction on 50 vertically oxidized structural configurations for each acene crystal. The results of our analysis are illustrated in Fig. 3(a) and (b) for naphthalene and pentacene, respectively. First, we note that, in both cases, average IPR values for vertically injected holes are close to those calculated for the ordered system, thus indicating that the disorder itself is not sufficient to properly localize the charge at a molecular level, and that calculations explicitly including the hole must be carried out to accurately describe the system. For the smaller acene, a vast majority of the investigated structural configurations feature IPR values closer to that corresponding to full localization on a single molecule. In fact, inspection of the lowest unoccupied molecular orbital (LUMO) indicates the formation of a molecular polaron with only a very minor degree of delocalization towards the first neighbours, cf. Fig. 3(c). At variance with this, in pentacene, two groups of MD configurations can be distinguished, based on the calculated IPR values: one bearing IPR values close to that of the semilocalized state close to the delocalized VB edge and another featuring polaronic localization, in a fashion similar to that observed for naphthalene, cf. Fig. 3(d). Overall, our analysis supports a physical picture in which, within the present model, room-temperature thermal disorder favours polaron formation even in the case of pentacene, for which this phenomenon appeared to be energetically unfavourable in the ordered crystal.⁴⁰

Fig. 3 Inverse participation ratio of the hole wave function as a function of the respective Kohn–Sham energy vs. the VBM of the neutral bulk in (a) naphthalene and (b) pentacene for 50 MD structural configurations. Vertical dashed lines represent the IPR associated with fully delocalized and localized charge states (cf. main text). Isodensity representation of the LUMO for a representative structural configuration in which the hole is localized almost entirely on a single (c) naphthalene and (d) pentacene molecule.

We note that literature reports on the spatial extents of large polarons in pentacene (e.g., ref. 101) arise from ESR analyses that, while compelling, depend on experimental conditions and on the model used to convert ESR lineshapes into a real-space extent. Subsequent ESR studies have demonstrated that carriers are localized at sufficiently low temperature and that motional-narrowing effects at higher temperature can mimic signatures of extended wavefunctions,^102–104 so an ESR-derived “molecular count” must be interpreted with care. Likewise, some simulation protocols that predict delocalization (for example FOB-SH and related nonadiabatic schemes) depends on approximations and some parameters,^17,105 such as classical nuclear dynamics (in this work, we also use classical MD, but we mitigated its potential limitations by using QM-derived force fields and refining the results with DFT calculations on top of MD snapshots), surface-hopping decoherence models, fragment basis and finite-size handling - that must be controlled to avoid an approximate estimate of delocalization degree, while in our case we do not have have any adjustable parameter. Taken together, these points motivate our qualitative comparison (single-molecule localization in naphtalene vs. mixed/semi-delocalized in pentacene) with experiment while acknowledging that both experimental and theoretical delocalization estimates depend on method and conditions.

Next, we discuss the polaron energy levels for naphthalene and pentacene at room temperature, cf. Fig. 4. These are calculated using the grand-canonical formulation of defects in crystalline semiconductors.^106,107 Within this theory, the formation energy of a defect X with charge q in the relaxed nuclear coordinates R_q is defined as follows:

where G[bulk] is the free energy of the pristine bulk system, G[X(R_q,q)] the free energy of the supercell bearing the defect with charge q, ε_v the VBM of the neutral bulk system, μ_i the chemical potential of the species i added/subtracted from the supercell n_i times, μ_e the chemical potential of the electron, and E_corr[X(R_q,q)] is an electrostatic finite-size correction accounting for the spurious interactions due to the presence of a charged periodic supercell (vide infra). Within this theory, the adiabatic charge transition level of a defect is defined as the electron chemical potential at which the formation free energy of the defect in charge state q is equal to that in the charge state q′.^106,107 i.e. :

Fig. 4 Schematic representation of the adiabatic hole polaron level (μ_ad(+/0)) and of the vertical reduction energy of the acene cation (μ_red(+/0)) in crystalline naphthalene and pentacene at room temperature (300 K), as calculated either for the ordered crystallographic structure or including disorder from the average value achieved from sampling of MD structural configurations. For each organic semiconductor, energies are referred to that of the VBM of the neutral ordered system.

Therefore, for a pristine acene crystal (i.e. without defects), the adiabatic energy level associated with the (+/0) charge transition level, i.e. the hole polaron level, reads as follows:¹⁰⁸

μ_ad(+/0) = −G[pol(R₊₁,+1)] + G[bulk] − E_corr[pol(R₊₁,+1)] − ε_v (6)

where G⁺[pol(R₊₁, +1)] is the free energy of the supercell bearing the localized hole polaron and E_corr[pol(R₊₁, +1)] the respective finite-size correction.

In DFT-based studies of defects/polarons, the free energies appearing in eqn (6) are usually approximated with the total energies of the relaxed defective and pristine systems at 0 K.¹⁰⁶ In contrast, when evaluating redox levels in aqueous solution or at the solid-liquid interface, molecular dynamics simulations combined with the thermodynamic integration methods are deployed to calculate the free-energy difference.^109–113 This technique connects the free-energy difference appearing in Eqn 6 with vertical energy differences.¹¹⁴ In particular, within the Marcus approximation,¹¹⁵

where 〈ΔE_ox〉_R0 is the energy difference associated with vertical oxidation of the neutral bulk system while 〈ΔE_ver〉_R+1 indicates the vertical reduction of the system containing the hole polaron, cf. Computational details section for details of the calculations. We note that, in our case, 〈ΔE_ver〉_R0 refers to vertical removal of an electron from a semi-localized band edge state, thus corresponding to sampling of the valence band edge of the system. Instead, 〈ΔE_ver〉_R+1 is associated with the vertical oxidation level of the localized polaron, which is defined as:

μ_ad(+/0) = 〈ΔE_ver〉_R+1 − ε_v + E_corr[pol(R₊₁,0)] − E_corr[pol(R₊₁,+1)] (8)

where E_corr[pol(R₊₁,0)] is the electrostatic finite-size correction needed to account for the effect of the ionic polarization charge, remaining in the vertically neutral supercell, upon removal of the hole from the system bearing the polaron.¹¹⁶

The correction term E_corr[pol(R₊₁, +1)] is calculated using the Freysoldt-Neugebauer-Van de Walle (FNV) scheme^106,107 commonly adopted for defects in crystalline semiconductors. This correction amounts to 0.21 and 0.08 eV for naphthalene and pentacene, respectively, considering the measured static dielectric constants of the two materials.^117,118 At variance with this, E_corr[pol(R₊₁,0)] is calculated using the formalism developed in ref. 119 and it is estimated to be 0.11 and 0.03 eV for naphthalene and pentacene, respectively.

When compared with its 0 K counterpart,⁴⁰ we observe for naphthalene that both the vertical reduction and the adiabatic charge transition levels are essentially shifted towards higher energies (approx. 0.18 eV) when including thermal disorder, in line with the shift of VBM induced by the broadening of the band at room temperature, as previously observed for inorganic semiconductors.³⁰ Therefore, when investigating a system in which the hole was already well localized in the gap (μ_ad = 0.36 eV and μ_red = 0.55 eV), the effect of thermal motion on the average level is simply a consequence of that experienced by the respective band edge, cf. SI for a comparison of the electronic DOS with and without thermal disorder. Then, we consider the levels for pentacene, cf. right panel in Fig. 4. For the perfectly ordered crystal, no significant charge localization has been envisaged (μ_ad = 0.02 eV, see also the DOS in SI). Furthermore, thermal renormalization of the VBM was found to be modest (0.03 eV).⁴⁰ We here calculate an average value of 0.09 (0.19) eV for the adiabatic (vertical) charge transition level, in accord with the envisaged molecular confinement of the hole (only structural configurations with high values of IPR have been considered in the average).

From the difference between vertical and adiabatic level, we can define the total reorganization energy associated with the localization:

λ = μ_red(+/0) − μ_ad(+/0). (9)

This quantity, in turn, can be split into the individual contributions of (i) the molecular reorganization λ_mol, i.e. energy gain upon relaxation of the molecule bearing the charge and (ii) the collective response of the surrounding environment, λ_env:

λ = λ_mol − λ_env. (10)

By computing λ_mol from gas-phase calculations, we can estimate the contribution of the medium in stabilizing the extra charge.

Values listed in Table 1 indicate that, while for both materials molecular reorganization represents the largest fraction of λ, as we calculate λ_mol = 0.11 (0.06) eV for naphthalene (pentacene), a significant contribution to stabilization is also provided by the environment, with λ_env = 0.08 (0.04) eV. The size of these quantities is relatively high when considering that acene crystals are apolar materials: they possess low static and high-frequency dielectric constants, and hence are usually expected to display a vanishingly small λ_env in continuum models.^120,121

Table 1 Calculated values of λ, λ_mol, and λ_env for hole polaron at room temperature, cf. main text for definitions. All values are given in eV

	λ	λmol	λenv
Naphthalene	0.19	0.11	0.08
Pentacene	0.10	0.06	0.04

---

Our values are somewhat larger than those evaluated in early studies based on classical molecular mechanics^42,43 even if the general trend is similar.⁴⁰ In stark contrast, recently a λ_env = 0.170 eV has been reported for pentacene, again based on classical molecular dynamics simulations,¹²² but employing a different methodology with respect to ref. 42 and 43, to account for statistical contributions. In this regard, we note that our force field-based molecular dynamics simulations do not include an explicitly polarizable force field. As such, the full electronic polarization of the crystal lattice in response to charge injection is not modeled during the dynamical phase. However, the re-evaluation of the total energies of MD structural configurations with DFT calculations including both a Koopmans' compliant fraction of Fock exchange and non-local electron correlations, should ensure an improved description of the energetics and, therefore, more accurate results. Moreover, for the apolar crystals with low dielectric constants, long-range polarization effects are expected to be limited. In such materials, the dominant contribution to environmental stabilization arises from structural reorganization of neighboring molecules, such as compression and reorientation. These features are explicitly captured in our MD trajectories (vide infra), see also Section S2.1 in the SI for a more detailed discussion, and our results are consistent with the pioneering work of Matyushov, who, employing a nonlocal response function theory, predicted non-negligible solvent reorganization energies for nonpolar media, and found correlation with the local rearrangement of neighbouring molecules towards the charged species, due to the strong electric field experienced by them.^44,45,123

Therefore, we conclude our analysis of the electronic structure of acene crystals in presence of thermal disorder, by surveying the structural features which can be related with λ_env. To this end, we evaluate the relative positioning (as achieved from MD simulations) of the reference molecule (either neutral or positively charged) with respect to the nearest neighbours. Since acene molecules are almost planar, the most important displacements can be described in terms of (i) distances between the centers of mass (r_CM) and (ii) Tait-Bryan angles (indicated as pitch, roll and yaw in the following, see Fig. 5), which are a standard notation to describe the orientation of a rigid body with respect to another plane, being used e.g. to describe the geometry of consecutive base pairs in DNA oligomers.^124,125

Fig. 5 Molecular orientation and nearest-neighbour symmetry in naphthalene and pentacene. Left: Schematic representation of Tait-Bryan Angles in a cartesian frame. Right: Representation of first-neighbours in naphthalene; each colour identifies symmetry-equivalent neighbours. For pentacene, the numbering is the same.

First of all, we note that the central molecule in the simulation cell has six nearest neighbours, which can be divided in three different couples depending on the symmetry of the system, as clearly illustrated in Fig. 5 (right) and in Fig. 1. In a fashion similar to other organic semiconductors, the associated crystal structures feature a high-mobility plane, whereas the mobility perpendicular to this plane is 1–2 orders of magnitude smaller.^18,126,127 Therefore, we here focus on the 2D-transport within the high-mobility plane, for which there can be up to three nearest-neighbor transfer integrals. Thus, allowing for some of the transfer integrals to be zero, or for pairs of parameters to be identical (which is the case of naphthalene, vide infra), such a lattice can describe the vast majority of high-mobility organic molecular semiconductors.

Thus, in the following, we will discuss the geometrical feature for each of the three symmetry-independent couples in the high-mobility plane (vide infra).

We start by focusing on the distance between the centers of mass. When comparing MD with either a charged or a neutral molecule in the center of the box, we notice that the distribution of distances denotes an appreciable difference in “solute–solvent” interactions when a cation is present, in particular for the stacked molecules (7–11) see Fig. 6, top panel and bottom panel. In fact, for both materials the overall r_CM distribution is moved towards lower distances. In particular, for naphthalene, we observe a shift of ≈0.1 Å in the peak and we also find the distribution to be noticeably broadened; instead, the curve for charged pentacene is found to feature two peaks, one closer to that of the neutral system and another at distances shorter by 0.14 Å. Such a trend is evident also when considering the other couples of interacting molecules (cf. SI), even if with minor displacements, due to the weaker intermolecular interactions in play with respect to the stacking couple.

Fig. 6 Distributions of r_CM (Å) between the neutral and charged reference molecule and its stacked neighbours for: naphthalene (a) and pentacene (b) crystals at 300 K. The curves are normalized so that the integral over all possible values is 1. The blue vertical dashed line indicates the value in the perfect crystal.

Most interestingly, when comparing the two acenes, we immediately notice that the distribution of distances and also of pitch, roll and yaw angles, are always significantly broader for naphthalene than for pentacene (see SI, Fig. S7–S14). This different behaviour can be rationalized by looking at the energy potential curves for movement around the equilibrium position (Fig. S2) along the stacking and angular directions described above: one immediately notices that the curves for pentacene show consistently lower energy minima and a much steeper profile, indicating that any fluctuation around equilibrium position is more hindered for pentacene than for naphthalene. This is consistent with the higher degree of localization observed for naphthalene and with the estimated values of reorganization energies. It also suggests that the effect of disorder on charge mobility will be different for the two systems. Overall, our analysis confirms that local rearrangements of the first neighbours surrounding the molecular polaron are at root of the non-negligible reorganization energies in an apolar environment, such as an acene crystal, in a fashion that could be considered similar to the electrostriction mechanism proposed within nonlocal response theory.^44,45,123

Fluctuations of electronic couplings

In addition to the energies, also the fluctuation in time of the electronic couplings plays an important role in organic semiconductors.

The electronic coupling element J_ij for charge transfer between two states i and j can be evaluated as

J_ij = 〈ϕ_i|F̂|ϕ_j〉 (11)

where ϕ_i and ϕ_j are the HOMO orbitals of the isolated monomers i and j, respectively, and F̂ is the Fock operator of the dimer.^{16,88,128,129}

We here evaluate the variance of the electronic coupling σ_ij²= 〈(J_ij − 〈J_ij〉)²〉, via the relation³⁶:

where S_ij,M (non-local electron-phonon coupling) is a measure of how the electronic coupling between molecular pair i and j is modulated by a displacement along a phonon mode M. This quantity can be evaluated by combining the Cartesian gradient of the transfer integral and the phonon modes of the system, as discussed in ref. 16, 80 and 81. More details about this procedure are reported in the Computational Details. Moreover, since the computation of crystal phonons is an extremely time-consuming process,^88,130 we have here lowered the computational burden by computing them resorting to the self-consistent charge density functional tight-binding (SCC-DFTB) method,⁸³ a non-orthogonal tight-binding method based on a second-order expansion of the DFT total-energy expression, which leads to computations about three-four orders of magnitude faster than DFT (see Computational details section).

Electronic couplings and their fluctuations, reported in Table 2, are in good agreement with previous reports.^{1,15,18,33,37} Inspection of that Table and Fig. 1 shows that, while naphthalene has two paths that are equivalent by symmetry, pentacene has three symmetry-nonequivalent paths. Moreover, |σ|/J is on average higher for naphthalene, with the exception of path C in pentacene. Since it is expected that (i) when |σ|/J is very large, the material should display a low mobility,^1,13,16 and (ii) the presence of two paths with lower |σ|/J (as in pentacene) ensures the possibility of circumventing a defect in the material, or of transferring the charge through percolation paths¹³¹ (e.g. in transistor or in bulk heterojunction solar cells), we expect that the disorder will have a higher impact in lowering the mobility for naphthalene than for pentacene, as we will discuss in the next section.

Table 2 Electronic coupling (J, absolute value), its fluctuation (σ), and their ratio, evaluated at T = 298 K for naphthalene and pentacene

	Path	J (meV)	σJ (meV)	σJ/J
Nap	A	45.4	25.4	0.56
	B	14.7	9.8	0.67
	B	14.7	9.8	0.67
Pn	A	40.2	19.6	0.49
	B	87.5	30.8	0.35
	C	44.6	34.3	0.77

---

Evaluation of charge mobility

It has been longly debated in the literature whether thermal fluctuations of the electronic coupling (nonlocal electron–phonon coupling, Peierls term),^1,12 or thermal site energy fluctuations (local electron–phonon couplings, Holstein term)^29,132 are the major source of polaron localization and mobility decay in organic crystals. To gain further insight on this issue, we calculate charge mobility for (i) a completely ordered crystal structure, or (ii) for a disordered system taking into account the spread in (a) electronic coupling and (b) polaron binding energies, modeled as a Gaussian distribution.

To this end, we employ Einstein–Smoluchowski's relation:

where:and n is the dimensionality of the system and 〈(r(t) − r(0))²〉 is the mean-square displacement.

Eqn (14) can be simplified for a perfect crystal structure.^15,27,133 However, we are here interested in the effect of energetic and coupling fluctuations and, therefore, we evaluate D by explicitly carrying out KMC simulations of charge transfer in a 20 × 20 × 1 supercell (since charge transport occurs mainly in the ab plane, we use a minimal replica along c direction). The influence of disorder on site energies and electronic couplings is accounted for by setting, at the start of each trajectory, the values of each site energy and each electronic coupling, randomly extracted from a Gaussian distribution. More details are given in the computational details section. The KMC trajectories are stopped after diffusive regime is reached, and we compute the diffusion coefficient D of the charge as in eqn (14), averaging over 20 000 independent KMC trajectories.

Charge transfer rates used in the KMC are evaluated on-the-fly resorting to Fermi's Golden Rule (FGR),

where F(ΔE_ij,T) is the Franck–Condon weighted density of states (FCWD) at the ΔE_ij of the transition and at a given temperature T, which is calculated from the geometry and frequencies of the neutral and charged states.^134–137 We chose FGR to evaluate the rates because, at variance with the widely used classical Marcus theory, the former allows taking into account tunneling effects and frequency changes between initial and final electronic states, including the whole sets of normal modes in the computation,^27,134–137 leading to accurate predictions of both radiative and non-radiative electronic transition rates.^{15,27,136,138}

We remark that extensions to the standard hopping KMC described above have been proposed in the literature. One of them resides on explicitly computing fluctuations of the electronic coupling of numerous dimers, extracted from MD fluctuations, as it has been done e.g. by Elstner group¹³⁹ and by some of us.⁹⁷ The drawback of this approach resides in the need for a comprehensive sampling of several molecule pairs, entailing a significant computational burden, which however could be alleviated by resorting to Machine Learning computations of couplings.⁹⁷ More advanced implementations of KMC have been recently introduced to deal with delocalized charge carriers typical of organic materials.^140,141 However, it has been shown that by using FGR the predictions are in line with the ones obtained resorting to a partially delocalized framework.⁹⁷ Thus, we will resort to the standard KMC approach with FGR rates, described in the previous paragraph.

The results in Table 3 show that neglecting the disorder leads to predictions strongly overestimated with respect to experimental data. As concerns the latter ones, a mobility around 1.0 cm² V⁻¹ s⁻¹ has been reported for naphthalene,^34,35,142 while μ ranging in the interval 0.6–2.3 cm² V⁻¹ s⁻¹ have been experimentally observed^38,143–146 for the pentacene polymorph^18,147,148 used in our analysis.¹⁴⁷ In contrast, when taking into account both the fluctuations of the energy and of the electronic couplings, we predict mobilities of 0.40 and 1.60 cm² V⁻¹ s⁻¹ for naphthalene and pentacene, respectively, in close agreement with experimental measurements.

Table 3 Calculated mobility (cm² V⁻¹ s⁻¹) neglecting (μ⁰) or taking into account the energetic fluctuations (μ^E), the fluctuations of electronic couplings (μ^J) or both (μ^D), for naphthalene and pentacene. Experimental mobilities (μ^exp)^{34,35,38,142–146} are also reported for comparison

	μ^0	μ^J	μ^E	μ^D	μ^exp
Nap	12.5	13.8	0.23	0.40	1.0
Pn	18.5	20.2	1.56	1.60	0.6–2.3

---

Importantly, we underline that not only the agreement with the experiment is retrieved when accounting for disorder, but also the difference in mobilities between naphthalene and pentacene is drastically influenced, as the energy levels of the former are much more affected. This indicates that neglecting the effect of disorder may lead to misleading results when carrying out comparative analyses of different OSs.

Finally, in order to disentangle the effect of the different factors on the final mobility, we have performed simulations in which either only energetic fluctuations or only electronic couplings fluctuations are taken into account. The results, reported in Table 3, show that the lowering in the mobility with respect to the perfect crystal structure is mainly related to the presence of site energy fluctuations, in line with what has been proposed in some works in the past.^29,132 Conversely, if only fluctuations of electronic couplings are taken into account, a slight increase in the mobility is predicted. Interestingly, also other models, e.g. transient localization theory^13,41 have led to similar results in the past, i.e. high fluctuations σ of the coupling value can have a beneficial impact on the charge transport properties of a molecule.¹⁵

The different impact of the two types of disorder on the mobility can be easily rationalized in light of the physical assumptions at the basis of the FGR/KMC model: indeed, a positive fluctuation in the coupling (i.e. an increase in the coupling), will lead to higher rates and thus increased mobility. On the other hand, energetic disorder can only lead to a strong decrease in the mobilities, since the FCWD for both the molecules (see Fig. S15 in the SI) is strongly peaked around E = 0, to the point that any fluctuation (either positive or negative) around this value will lead to a decrease in the FCWD and thus (eqn 3) in the rates determining the charge mobility. However, it must be kept in mind that the opposing roles of diagonal (energetic) and off-diagonal (coupling) disorder have a clear physical origin and are not simply artefacts of neglecting coherence. Diagonal disorder generates energetic traps and Anderson-type localization that reduce rates in both incoherent (hopping) and coherent pictures; conversely, fluctuating transfer integrals can raise the average transition rate, that further enhances mobility.

The broadening of energy levels and electronic couplings in OSs, due to their soft nature has been largely discussed in literature^1,12,16 along with their impact on charge mobilities. While some studies highlighted the relationship between fluctuations of the electronic coupling and the calculated rates (e.g. ref. 1 and 12), other works have also discussed the relevance of disorder on the site energy (e.g. ref. 13, 29 and 132). However, these studies mainly focused on the struggle between localization and delocalization of the charge carrier, with thermal disorder favouring polaron formation and, consequently, slowing down transport. While this is qualitatively consistent with the physical picture ensuing from our simulations, i.e. disorder promoting and stabilizing localization, we here show that the spread of polaronic binding energies - that is, the distribution of charge stabilization energies across thermally disordered configurations - generated by intramolecular and intermolecular motion (local rearrangement induced by the positively charged molecule on the surrounding molecules), needs to be accurately accounted for when describing charge mobility: if only the average polaron binding energy were considered, such inhomogeneities would be neglected, resulting in an overly idealized transport model that would likely overestimate mobilities (as in our case). Therefore, including the full distribution of binding energies is crucial to reproduce the experimental trends and to correctly capture the disorder-driven mechanistic origin of transport suppression. In this respect, it must be remarked that, since this effect is material-dependent, it cannot be ignored in screening procedures.

Conclusions

In this work, we combined classical MD simulations based on quantum-mechanically derived force fields and DFT calculations based on a Koopmans' compliant formulation of hybrid functionals to investigate the electronic structure of a hole polaron in naphthalene and pentacene at room temperature. Our analysis demonstrates that thermal disorder shifts the polaronic hole energy levels towards higher energy, thus stabilizing the localized charge. In particular, room-temperature disorder allows for polaron formation even in the case of pentacene, for which no localization is predicted when considering the ordered crystal.

The overall stabilization of the polaron, while mainly stemming from molecular reorganization of a single acene, is found to feature also a significant contribution from the environment, in line with prediction based on non-local response theory.⁴⁴ Furthermore, by inspection of the “micro-solvation” of the molecularly localized polaron, we are able to identify the structural features associated with the calculated reorganization energies, namely a contraction of the first neighbours towards the positive charge, i.e. electrostriction, and noticeable distortions in the inter-molecular interactions, with respect to the neutral system.

We emphasize that our results are obtained using a state-of-the-art modeling framework, which however relies on two key assumptions. First, we assume that classical MD, based on quantum-mechanically derived force fields, can accurately describe the distribution of static defects in the crystal. We consider this assumption well justified, since in previous work on related systems we have found excellent agreement between quantum and classical (with QM-derived force fields) MD.⁴⁰ Additionally, we have refined classical MD configurations through DFT optimization, further mitigating potential inaccuracies. In any case, we are working to compare the classical MD reported here and a full MD performed at the hybrid DFT level in a forthcoming work.

Second, we assume that charge transport can be reliably described using kinetic Monte Carlo simulations based on a localized basis. This assumption is consistent with the results of DFT calculations on MD snapshots, which point toward a small polaron model. In addition, the predicted FGR hole hopping rates are not unphysically high, but range from tens of fs to a few ps, so that adopting full quantum dynamic calculations,¹²⁹ which are able to account for coherent effects, appear to be not necessary here.

Noteworthy, our FGR approach goes well beyond classical Marcus theory, including tunneling effects, which are known to improve the realism of charge transport predictions.^27,28

By combining classical MD sampling, DFT optimization, and KMC simulations over 20 000 trajectories, our methodology captures both site-energy and electronic coupling disorder, yielding charge mobilities in excellent agreement with experimental measurements for both naphthalene and pentacene. Further analyses show that fluctuations of site energies have a significant impact on charge mobilities, lowering them by ca. one order of magnitude, whereas electronic coupling fluctuations appear to slightly facilitate charge transport, in line with what has been suggested in previous works.^29,132,149

Overall, this work demonstrates that, within the validity range of our methodology - namely the absence of ultrafast transitions ensuring the applicability of FGR, and the use of QM-derived force fields to reliably describe the system's time evolution - our modeling approach provides a physically grounded and computationally efficient protocol for understanding polaron formation and charge transport in thermally disordered organic semiconductors.

Author contributions

A. Landi, F.A. and A.P. contributed equally to the conceptualization of the study. A. Landi led the formal analysis and software development, and shared responsibilities in data curation, investigation, methodology, supervision, and writing of the original draft and subsequent revisions. F. A. contributed to formal analysis, investigation, methodology, and co-authored both the original draft and revisions. A. Leo and D. P. participated in the investigation and review of the manuscript, with D. P. also contributing to software development. G. P. was involved in formal analysis, investigation, software development, and co-authorship of the original draft and revisions. A. P. contributed to the conceptualization, project administration, supervision, and review of the manuscript.

Conflicts of interest

There are no conflicts to declare.

Data availability

The data supporting this article have been included as part of the Supplementary information (SI). Supplementary Information: additional computational details concerning QMD–FFs parameterisation and validation. Density of states for neutral and positively charged supercells. Details about analysis of distances and angles distribution for interacting molecules from MD simulations. Discussion about different theoretical models for charge transport. FCWD for naphthalene and pentacene. Details for KMC simulation protocol. The force fields used in this work have been made available to the community in a public GitHub repository (https://github.com/AlessandroUniSa1/acenes_paper2025). See DOI: https://doi.org/10.1039/d5tc03115j.

Acknowledgements

A.Landi gratefully acknowledges PRIN 2022 grant 2022XSC9P5 (INTESA-SOLE) and PNRR MUR Project CN00000013-ICSC (INNOVATOR) by Ministero of Università e Ricerca (MUR) for funding. A.Landi acknowledges the CINECA award HP10CVHRXV (LIMES) under the ISCRA initiative, for the availability of high performance computing resources and support. F.A. thankfully acknowledges REACT-EU DM 1062/2021 for funding and the CINECA projects AID-Loc, Def-OCS, and Reo-Solv for computational resources. A.P. gratefully acknowledges PRIN PNRR 2022 grant 2022WXPMB. D.P. gratefully acknowledges funding from the Italian Ministry of University and Research (MUR, PRIN grant 2022CXHY3A) and from Spoke 7 of ICSC Centro Nazionale di Ricerca in HPC, Big Data e Quantum Computing within PNRR.

References

1. S. Fratini, S. Ciuchi, D. Mayou, G. T. De Laissardière and A. Troisi, Nat. Mater., 2017, 16, 998–1002.
2. E. K. Solak and E. Irmak, RSC Adv., 2023, 13, 12244–12269.
3. T. F. O'Connor, A. V. Zaretski, S. Savagatrup, A. D. Printz, C. D. Wilkes, M. I. Diaz, E. J. Sawyer and D. J. Lipomi, Sol. Energy Mater. Sol. Cells, 2016, 144, 438–444.
4. T. Kubo, R. Häusermann, J. Tsurumi, J. Soeda, Y. Okada, Y. Yamashita, N. Akamatsu, A. Shishido, C. Mitsui, T. Okamoto, S. Yanagisawa, H. Matsui and J. Takeya, Nat. Commun., 2016, 7, 11156.
5. A. Landi, G. Ricci, Y. Olivier, A. Capobianco and A. Peluso, J. Phys. Chem. Lett., 2024, 15, 11042–11050.
6. K. Liu, B. Ouyang, X. Guo, Y. Guo and Y. Liu, npj Flex. Electron., 2022, 6, 1.
7. A. Landi, M. Reisjalali, J. D. Elliott, M. Matta, P. Carbone and A. Troisi, J. Mater. Chem. C, 2023, 11, 8062–8073.
8. B. D. Paulsen, K. Tybrandt, E. Stavrinidou and J. Rivnay, Nat. Mater., 2020, 19, 13–26.
9. O. Ostroverkhova, Chem. Rev., 2016, 116, 13279–13412.
10. K. Yoshikawa, H. Kawasaki, W. Yoshida, T. Irie, K. Konishi, K. Nakano, T. Uto, D. Adachi, M. Kanematsu, H. Uzu and K. Yamamoto, Nat. Ener., 2017, 2, 17032.
11. K. Liu, B. Ouyang, X. Guo, Y. Guo and Y. Liu, npj Flex. Electron., 2022, 6, 1.
12. A. Troisi and G. Orlandi, Phys. Rev. Lett., 2006, 96, 086601.
13. S. Fratini, D. Mayou and S. Ciuchi, Adv. Funct. Mater., 2016, 26, 2292–2315.
14. T. Nematiaram and A. Troisi, J. Chem. Phys, 2020, 152, 190902.
15. A. Landi, J. Phys. Chem. C, 2019, 123, 18804–18812.
16. A. Landi and A. Troisi, J. Phys. Chem. C, 2018, 122, 18336–18345.
17. S. Giannini, A. Carof, M. Ellis, H. Yang, O. G. Ziogos, S. Ghosh and J. Blumberger, Nat. Commun., 2019, 10, 3843.
18. A. Landi, A. Troisi and A. Peluso, J. Mater. Chem. C, 2019, 7, 9665–9670.
19. S. Fratini, M. Nikolka, A. Salleo, G. Schweicher and H. Sirringhaus, Nat. Mater., 2020, 19, 491–502.
20. A. Troisi, Adv. Mater., 2007, 19, 2000–2004.
21. F. Ortmann, F. Bechstedt and K. Hannewald, Phys. Status Solidi B, 2011, 248, 511–525.
22. J. Ren, N. Vukmirović and L.-W. Wang, Phys. Rev. B: Condens. Matter Mater. Phys., 2013, 87, 205117.
23. A. Heck, J. J. Kranz, T. Kubar and M. Elstner, J. Chem. Theory Comput., 2015, 11, 5068–5082.
24. W. Xie, D. Holub, T. Kubar and M. Elstner, J. Chem. Theory Comput., 2020, 16, 2071–2084.
25. S. Fratini, S. Ciuchi, D. Mayou, G. T. De Laissardière and A. Troisi, Nat. Mater., 2017, 16, 998–1002.
26. Y. C. Cheng, R. J. Silbey, D. A. da Silva Filho, J. P. Calbert, J. Cornil and J. L. Brédas, J. Chem. Phys, 2003, 118, 3764–3774.
27. L. Wang, G. Nan, X. Yang, Q. Peng, Q. Li and Z. Shuai, Chem. Soc. Rev., 2010, 39, 423–434.
28. G. Nan, X. Yang, L. Wang, Z. Shuai and Y. Zhao, Phys. Rev. B: Condens. Matter Mater. Phys., 2009, 79, 115203.
29. S. Giannini, A. Carof and J. Blumberger, J. Phys. Chem. Lett., 2018, 9, 3116–3123.
30. J. Wiktor, F. Ambrosio and A. Pasquarello, ACS Energy Lett., 2018, 3, 1693–1697.
31. S. Giannini, P. M. Martinez, A. Semmeq, J. P. Galvez, A. Piras, A. Landi, D. Padula, J. G. Vilhena, J. Cerezo and G. Prampolini, J. Chem. Theory Comput., 2025, 21, 3156–3175.
32. J. Cerezo, G. Prampolini and I. Cacelli, Theor. Chem. Acc., 2018, 137, 80.
33. R. S. Sánchez-Carrera, P. Paramonov, G. M. Day, V. Coropceanu and J.-L. Brédas, J. Am. Chem. Soc., 2010, 132, 14437–14446.
34. N.-E. Lee, J.-J. Zhou, L. A. Agapito and M. Bernardi, Phys. Rev. B, 2018, 97, 115203.
35. W. Warta and N. Karl, Phys. Rev. B: Condens. Matter Mater. Phys., 1985, 32, 1172–1182.
36. V. Coropceanu, R. S. Sánchez-Carrera, P. Paramonov, G. M. Day and J.-L. Brédas, J. Phys. Chem. C, 2009, 113, 4679–4686.
37. Y. Yi, V. Coropceanu and J.-L. Brédas, J. Chem. Phys, 2012, 137, 164303.
38. Y.-Y. Lin, D. I. Gundlach, S. F. Nelson and T. N. Jackson, IEEE Trans. Electron Devices, 1997, 44, 1325–1331.
39. T. Sekitani, Y. Kato, S. Iba, H. Shinaoka, T. Someya, T. Sakurai and S. Takagi, Appl. Phys. Lett, 2005, 86, 073511.
40. F. Ambrosio, J. Wiktor, A. Landi and A. Peluso, J. Phys. Chem. Lett., 2023, 14, 3343–3351.
41. S. Ciuchi, R. C. Hatch, H. Höchst, C. Faber, X. Blase and S. Fratini, Phys. Rev. Lett., 2012, 108, 256401.
42. D. P. McMahon and A. Troisi, J. Phys. Chem. Lett., 2010, 1, 941–946.
43. J. E. Norton and J.-L. Brédas, J. Am. Chem. Soc., 2008, 130, 12377–12384.
44. D. V. Matyushov, Phys. Chem. Chem. Phys., 2020, 22, 10653–10665.
45. D. V. Matyushov, Phys. Chem. Chem. Phys., 2023, 25, 7589–7610.
46. R. A. Holroyd and J. R. Miller, J. Phys. Chem. B, 2019, 123, 9206–9211.
47. G. Ricci, S. Canola, Y. Dai, D. Fazzi and F. Negri, Molecules, 2021, 26, 4119.
48. J. Lederer, W. Kaiser, A. Mattoni and A. Gagliardi, Adv. Theory Simul., 2019, 2, 1800136.
49. S. M. Gali, G. D’Avino, P. Aurel, G. Han, Y. Yi, T. A. Papadopoulos, V. Coropceanu, J.-L. Brédas, G. Hadziioannou, C. Zannoni and L. Muccioli, J. Chem. Phys., 2017, 147, 134904.
50. C. R. Groom, I. J. Bruno, M. P. Lightfoot and S. C. Ward, Acta Crystallogr., Sect. B: Struct. Sci., Cryst. Eng. Mater., 2016, 72, 171–179.
51. A. Landi and D. Padula, J. Mater. Chem. A, 2021, 9, 24849–24856.
52. W. L. Jorgensen, D. S. Maxwell and J. Tirado-Rives, J. Am. Chem. Soc., 1996, 118, 11225–11236.
53. C. I. Bayly, P. Cieplak, W. Cornell and P. A. Kollman, J. Phys. Chem., 1993, 97, 10269–10280.
54. M. J. Abraham, T. Murtola, R. Schulz, S. Páll, J. C. Smith, B. Hess and E. Lindahl, SoftwareX, 2015, 1–2, 19–25.
55. T. Darden, D. York and L. Pedersen, J. Chem. Phys, 1993, 98, 10089–10092.
56. B. Hess, H. Bekker, H. J. C. Berendsen and J. G. E. M. Fraaije, J. Comput. Chem., 1997, 18, 1463–1472.
57. H. J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren, A. DiNola and J. R. Haak, J. Chem. Phys, 1984, 81, 3684–3690.
58. M. Parrinello and A. Rahman, J. Appl. Phys, 1981, 52, 7182–7190.
59. J. P. Perdew, M. Ernzerhof and K. Burke, J. Chem. Phys, 1996, 105, 9982–9985.
60. C. Adamo and V. Barone, J. Chem. Phys, 1999, 110, 6158–6170.
61. G. Miceli, W. Chen, I. Reshetnyak and A. Pasquarello, Phys. Rev. B, 2018, 97, 121112.
62. T. Bischoff, I. Reshetnyak and A. Pasquarello, Phys. Rev. B, 2019, 99, 201114.
63. T. Bischoff, J. Wiktor, W. Chen and A. Pasquarello, Phys. Rev. Mater., 2019, 3, 123802.
64. T. Bischoff, I. Reshetnyak and A. Pasquarello, Phys. Rev. Research, 2021, 3, 023182.
65. J. P. Perdew and A. Zunger, Phys. Rev. B: Condens. Matter Mater. Phys., 1981, 23, 5048–5079.
66. Y. Zhang and W. Yang, J. Chem. Phys, 1998, 109, 2604–2608.
67. O. A. Vydrov and T. Van Voorhis, J. Chem. Phys, 2010, 133, 244103.
68. R. Sabatini, T. Gorni and S. de Gironcoli, Phys. Rev. B: Condens. Matter Mater. Phys., 2013, 87, 041108.
69. J. VandeVondele, M. Krack, F. Mohamed, M. Parrinello, T. Chassaing and J. Hutter, Comput. Phys. Commun., 2005, 167, 103–128.
70. C. Hartwigsen, S. Goedecker and J. Hutter, Phys. Rev. B: Condens. Matter Mater. Phys., 1998, 58, 3641–3662.
71. J. VandeVondele and J. Hutter, J. Chem. Phys, 2007, 127, 114105.
72. M. Guidon, F. Schiffmann, J. Hutter and J. VandeVondele, J. Chem. Phys., 2008, 128, 214104.
73. M. Guidon, J. Hutter and J. VandeVondele, J. Chem. Theory Comput., 2010, 6, 2348–2364.
74. A. Kubas, F. Gajdos, A. Heck, H. Oberhofer, M. Elstner and J. Blumberger, Phys. Chem. Chem. Phys., 2015, 17, 14342–14354.
75. T. P. Nguyen, J. H. Shim and J. Y. Lee, J. Phys. Chem. C, 2015, 119, 11301–11310.
76. D. Padula, A. Landi and G. Prampolini, Energy Adv., 2023, 2, 1215–1224.
77. S. Miertuš, E. Scrocco and J. Tomasi, Chem. Phys., 1981, 55, 117–129.
78. G. D'Avino, L. Muccioli, C. Zannoni, D. Beljonne and Z. G. Soos, J. Chem. Theory Comput., 2014, 10, 4959–4971.
79. 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.
80. X. Xie, A. Santana-Bonilla and A. Troisi, J. Chem. Theory Comput., 2018, 14, 3752–3762.
81. A. Landi, A. Peluso and A. Troisi, Adv. Mater., 2021, 33, 2008049.
82. B. Hourahine, B. Aradi, V. Blum, F. Bonafé, A. Buccheri, C. Camacho, C. Cevallos, M. Y. Deshaye, T. Dumitrică, A. Dominguez, S. Ehlert, M. Elstner, T. van der Heide, J. Hermann, S. Irle, J. J. Kranz, C. Köhler, T. Kowalczyk, T. Kubař, I. S. Lee, V. Lutsker, R. J. Maurer, S. K. Min, I. Mitchell, C. Negre, T. A. Niehaus, A. M. N. Niklasson, A. J. Page, A. Pecchia, G. Penazzi, M. P. Persson, J. Řezáč, C. G. Sánchez, M. Sternberg, M. Stöhr, F. Stuckenberg, A. Tkatchenko, V. W.-Z. Yu and T. Frauenheim, J. Chem. Phys., 2020, 152, 124101.
83. M. Elstner, D. Porezag, G. Jungnickel, J. Elsner, M. Haugk, T. Frauenheim, S. Suhai and G. Seifert, Phys. Rev. B: Condens. Matter Mater. Phys., 1998, 58, 7260–7268.
84. M. Gaus, A. Goez and M. Elstner, J. Chem. Theory Comput., 2013, 9, 338–354.
85. H. J. Monkhorst and J. D. Pack, Phys. Rev. B, 1976, 13, 5188–5192.
86. M. T. Ruggiero, S. Ciuchi, S. Fratini and G. D’Avino, J. Phys. Chem. C, 2019, 123, 15897–15907.
87. A. K. Rappe, C. J. Casewit, K. S. Colwell, W. A. Goddard and W. M. Skiff, J. Am. Chem. Soc., 1992, 114, 10024–10035.
88. T. Nematiaram, D. Padula, A. Landi and A. Troisi, Adv. Funct. Mater., 2020, 30, 2001906.
89. R. Borrelli and A. Peluso, MolFC: A program for Franck-Condon integrals calculation, Package available online at https://www.theochem.unisa.it.
90. R. Borrelli, M. Di Donato and A. Peluso, Biophys. J., 2005, 89, 830–841.
91. A. Landi, A. Landi, A. Leo and A. Peluso, J. Chem. Phys, 2024, 160, 174114.
92. A. Capobianco, J. Wiktor, A. Landi, F. Ambrosio and A. Peluso, Nano Lett., 2024, 24, 8335–8342.
93. A. Peluso, R. Borrelli and A. Capobianco, J. Phys. Chem. A, 2013, 117, 10985.
94. Y. Niu, W. Li, Q. Peng, H. Geng, Y. Yi, L. Wang, G. Nan, D. Wang and Z. Shuai, Mol. Phys., 2018, 116, 1078–1090.
95. A. Landi, A. Landi, A. Velardo and A. Peluso, ACS Appl. Energy Mater., 2022, 5, 10815–10824.
96. H. van Eersel, R. A. J. Janssen and M. Kemerink, Adv. Funct. Mater., 2012, 22, 2700–2708.
97. D. Padula, L. Barneschi, A. Peluso, T. Cinaglia and A. Landi, J. Mater. Chem. C, 2023, 11, 12297–12306.
98. P. Reiser, M. Konrad, A. Fediai, S. Léon, W. Wenzel and P. Friederich, J. Chem. Theory Comput., 2021, 17, 3750–3759.
99. V. Rühle, A. Lukyanov, F. May, M. Schrader, T. Vehoff, J. Kirkpatrick, B. Baumeier and D. Andrienko, J. Chem. Theory Comput., 2011, 7, 3335–3345.
100. K. Kaur and C. V. Singh, Energy Procedia, 2012, 29, 291–299.
101. K. Marumoto, S.-i Kuroda, T. Takenobu and Y. Iwasa, Phys. Rev. Lett., 2006, 97, 256603.
102. H. Matsui, A. S. Mishchenko and T. Hasegawa, Phys. Rev. Lett., 2010, 104, 056602.
103. H. Matsui, T. Hasegawa, Y. Tokura, M. Hiraoka and T. Yamada, Phys. Rev. Lett., 2008, 100, 126601.
104. A. S. Mishchenko, H. Matsui and T. Hasegawa, Phys. Rev. B: Condens. Matter Mater. Phys., 2012, 85, 085211.
105. H. Oberhofer, K. Reuter and J. Blumberger, Chem. Rev., 2017, 117, 10319–10357.
106. C. Freysoldt, J. Neugebauer and C. G. Van de Walle, Phys. Rev. Lett., 2009, 102, 016402.
107. H.-P. Komsa, T. T. Rantala and A. Pasquarello, Phys. Rev. B: Condens. Matter Mater. Phys., 2012, 86, 045112.
108. F. Ambrosio, J. Wiktor, F. De Angelis and A. Pasquarello, Energy Environ. Sci., 2018, 11, 101–105.
109. J. Cheng and M. Sprik, Phys. Chem. Chem. Phys., 2012, 14, 11245–11267.
110. F. Ambrosio, G. Miceli and A. Pasquarello, J. Chem. Phys, 2015, 143, 244508.
111. F. Ambrosio, J. Wiktor and A. Pasquarello, ACS Appl. Mater. Interfaces, 2018, 10, 10011–10021.
112. Z. Guo, F. Ambrosio and A. Pasquarello, ACS Catalysis, 2020, 10, 13186–13195.
113. G. Di Liberto, F. Maleki and G. Pacchioni, J. Phys. Chem. C, 2022, 126, 10216–10223.
114. D. Frenkel and B. Smit, Understanding Molecular Simulation: From Algorithms to Applications, Academic Press, 2002.
115. F. Ambrosio, A. Landi, A. Peluso and A. Capobianco, J. Chem. Theory Comput., 2024, 20, 9708–9719.
116. S. Falletta, J. Wiktor and A. Pasquarello, Phys. Rev. B, 2020, 102, 041115.
117. C. H. Kim, O. Yaghmazadeh, D. Tondelier, Y. B. Jeong, Y. Bonnassieux and G. Horowitz, J. Appl. Phys, 2011, 109, 083710.
118. K. Ishii, M. Kinoshita and H. Kuroda, Bull. Chem. Soc. Jpn., 1973, 46, 3385–3391.
119. S. Falletta and A. Pasquarello, Phys. Rev. B, 2023, 107, 205125.
120. R. A. Marcus, J. Chem. Phys, 1956, 24, 966–978.
121. R. A. Marcus, J. Chem. Phys, 1956, 24, 979–989.
122. C.-H. Yang, C.-I. Wang, Y.-S. Wang and C.-P. Hsu, J. Chem. Theory Comput., 2024, 20, 6981–6991.
123. G. Kodis, I. R. Gould and D. V. Matyushov, Phys. Rev. Res., 2021, 3, 013109.
124. X.-J. Lu and W. K. Olson, Nat. Protoc., 2008, 3, 1213–1227.
125. A. Capobianco, A. Landi and A. Peluso, Phys. Chem. Chem. Phys., 2017, 19, 13571–13578.
126. W. Deng, L. Sun, J. Huang, S. Chai, S. Wen and K. Han, Nature Protocols, 2015, 10, 632–642.
127. B. Blülle, A. Troisi, R. Häusermann and B. Batlogg, Phys. Rev. B, 2016, 93, 035205.
128. A. Troisi and G. Orlandi, Chem. Phys. Lett., 2001, 344, 509–518.
129. A. Landi, A. Capobianco and A. Peluso, J. Phys. Chem. Lett., 2020, 11, 7769–7775.
130. T. F. Harrelson, V. Dantanarayana, X. Xie, C. Koshnick, D. Nai, R. Fair, S. A. Nuñez, A. K. Thomas, T. L. Murrey, M. A. Hickner, J. K. Grey, J. E. Anthony, E. D. Gomez, A. Troisi, R. Faller and A. J. Moulé, Mater. Horiz., 2019, 6, 182–191.
131. J. S. Mehta and J. M. Mativetsky, ACS Appl. Energy Mater., 2018, 1, 5656–5662.
132. L. Wang and D. Beljonne, J. Phys. Chem. Lett., 2013, 4, 1888–1894.
133. W.-Q. Deng and W. A. Goddard, J. Phys. Chem. B, 2004, 108, 8614–8621.
134. R. Kubo and Y. Toyozawa, Prog. Theor. Phys., 1955, 13, 160–182.
135. M. Lax, J. Chem. Phys, 1952, 20, 1752–1760.
136. A. Landi, R. Borrelli, A. Capobianco, A. Velardo and A. Peluso, J. Chem. Theory Comput., 2018, 14, 1594–1601.
137. R. Borrelli and A. Peluso, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2013, 3, 542–559.
138. A. Velardo, R. Borrelli, A. Capobianco, A. Landi and A. Peluso, J. Phys. Chem. C, 2019, 123, 14173–14179.
139. A. D. Özdemir, S. Inanlou, F. Symalla, W. Xie, W. Wenzel and M. Elstner, J. Chem. Theory Comput., 2023, 19, 3849–3860.
140. D. Balzer, T. J. A. M. Smolders, D. Blyth, S. N. Hood and I. Kassal, Chem. Sci., 2021, 12, 2276–2285.
141. D. Balzer and I. Kassal, J. Phys. Chem. Lett., 2023, 2155–2162.
142. W. Warta, R. Stehle and N. Karl, Appl. Phys. A: Solids Surf., 1985, 36, 163–170.
143. J. Y. Lee, S. Roth and Y. W. Park, Appl. Phys. Lett, 2006, 88, 252106.
144. T. Uemura, M. Yamagishi, J. Soeda, Y. Takatsuki, Y. Okada, Y. Nakazawa and J. Takeya, Phys. Rev. B: Condens. Matter Mater. Phys., 2012, 85, 035313.
145. L. B. Roberson, J. Kowalik, L. M. Tolbert, C. Kloc, R. Zeis, X. Chi, R. Fleming and C. Wilkins, J. Am. Chem. Soc., 2005, 127, 3069–3075.
146. H. Yang, T. J. Shin, M.-M. Ling, K. Cho, C. Y. Ryu and Z. Bao, J. Am. Chem. Soc., 2005, 127, 11542–11543.
147. R. B. Campbell, J. M. Robertson and J. Trotter, Acta Cryst., 1962, 15, 289–290.
148. C. C. Mattheus, G. A. de Wijs, R. A. de Groot and T. T. M. Palstra, J. Am. Chem. Soc., 2003, 125, 6323–6330.
149. X. Bai, J. Qiu and L. Wang, J. Chem. Phys, 2018, 148, 104106.

---