Modeling the electronic structure of organic materials: A solid-state physicist's perspective

Modeling the electronic and optical properties of organic semiconductors remains a challenge for theory, despite the remarkable progress achieved in the last three decades. The complexity of these systems, including structural (dis)order and the still debated doping mechanisms, has been engaging theorists with different backgrounds. Regardless of the common interest across the various communities active in this field, these efforts have not led so far to a truly interdisciplinary research area. In the attempt to move further in this direction, we present our perspective as solid-state theorists for the study of molecular materials in different states of matter. Considering exemplary systems belonging to well-known families of oligo-acenes and -thiophenes, we provide a quantitative description of electronic properties and optical excitations obtained with state-of-the-art first-principles methods such as density-functional theory and many-body perturbation theory. Simulating the systems as gas-phase molecules, clusters, and periodic lattices, we are able to identify short- and long-range effects in their electronic structure. While the latter are usually dominant in organic crystals, the former play an important role, too, especially in the case of donor/acceptor complexes. Furthermore, we demonstrate the viability of implicit schemes to evaluate band gaps of molecules embedded in isotropic and even anisotropic environments, in quantitative agreement with experiments. In the context of doped organic semiconductors, we show how the crystalline packing enhances the favorable characteristics of these systems for opto-electronic applications. The counter-intuitive behavior predicted for their electronic and optical properties is deciphered with the aid of a tight-binding model, which represents a connection to the most common approaches to evaluate transport properties in these materials.


Introduction
Organic semiconductors are an established class of materials [1] that have received considerable attention in the last three decades by a vast community across physics, chemistry, materials science, and engineering [2]. Since the first successes in the synthesis of ordered films of organic molecules and their integration in nanotransistors [3], the interest in organic semiconductors as opto-electronic materials has exploded. Their structural flexibility, chemical tunability, and low production costs have attracted hundreds of scientists and technologists to the dream of plastic electronic devices [4,5]. Such a formidable experimental interest has driven the solid-state theorists' community to look systematically into the electronic and optical properties of these materials. Pioneering results obtained from first principles date back to the last decade of the 20th century [6] and to the early 2000s [7,8,9,10,11,12,13,14]. These works have been largely complemented by studies employing model Hamiltonians, especially for the investigation of complex phenomena involving electron-phonon coupling [15,16,17,18,19,20,21,22]. In parallel, the community active in quantum chemistry engaged itself in the investigation of organic materials by applying highly accurate methods such as coupled-cluster and configuration interaction [23,24,25,26], as well as efficient semiempirical schemes based on the Hartree-Fock approximation [27,28,29]. While both groups largely contributed to the interpretation of the experimental findings and to the progress of this research line, their efforts have remained mainly confined within their respective communities.
In the second decade of the 21st century, the hype for organic semiconductors was substantially reduced by the rise of low-dimensional semiconductors [30] and halide perovskites [31] as novel materials for opto-electronic applications. These emerging research areas absorbed most of the interest of the scientific community, including theorists. Specifically, the first-principle investigation of hybrid halide perovskites has appeared to be exceptionally challenging from the very beginning. The accurate description of their electronic structure demands the inclusion of spin-orbit coupling [32,33,34]. Moreover, the complex structure of these materials combined with the peculiar softness of their lattices requires non-standard methods [35,36,37,38,39]: for example, it was recently demonstrated that the lattice contribution to the screened Coulomb interaction should be explicitly incorporated in order to obtain accurate estimates of the exciton binding energies [40].
In such a context, molecular semiconductor modelling was devoted to addressing more specific questions related, among others, to doping mechanisms [41,42,43,44,45,46,47], polymorphism [48,49,50,51,52], and singlet fission [53,54,55,56,57,58,59]. Moreover, increasing attention was dedicated to disclosing the excitonic properties of organic materials [52,60,61,62], including quantifying resonance energies and clarifying the character of electron-hole pairs. These efforts, requiring an exceptionally high degree of specialization, enhanced the segregation between quantum/computational chemistry and solid-state theory. Both communities have relevantly contributed to address questions related to the role of the environment [63,64,46,65] and of electronvibrational couplings [66,67,68] in the opto-electronic activity of organic materials, as well as to benchmark the performance of advanced first-principles methodologies applied to this material class [69,70,61]. Yet, endeavors in physics and chemistry have mostly run in parallel with very rare intersections. Only by merging the efforts of different communities it will be possible to solve problems of increasing complexity leading to a substantial advancement in the fundamental and applied research on organic materials.
Adopting our perspective as solid-state physicists, we discuss herein the ability of density-functional theory and many-body perturbation theory, in their respective implementations for isolated and extended systems, to describe the electronic and optical properties of organic materials, ranging from isolated molecules to clusters, and crystalline solids. Taking anthracene as an example, we assess the role of short-and long-range effects in electronic structure and optical excitations. Given the numerical complexity of full-fledged simulations on organic crystals, we prove the viability of implicit schemes coupled to density-functional theory to account for electrostatic interactions exerted on molecules by both isotropic and anisotropic environments, obtaining band-gap values in excellent agreement with experiments. Finally, turning to donor/acceptor complexes, we discuss the importance of accounting for the lattice periodicity in order to describe the more favorable opto-electronic properties of cocrystals in comparison with their non-periodic counterparts. To interpret the counterintuitive behavior of these systems, we complement our first-principles results with a tight-binding model, thus providing a connection with established schemes to describe transport properties.
This review is organized as follows: In Section 2, we introduce the basic formalism of density-functional theory (Subsection 2.1) and many-body perturbation theory (Subsection 2.2). Next, we discuss the electronic and optical properties of anthracene, first described atomistically in different states of matter (Section 3) and then simulated as a single molecule embedded in various environments (Section 4). Finally, we examine the physics of a well-known donor/acceptor complex formed by a p-doped quaterthiophene oligomer, making use of a tight-binding model to shed light on the electronic and optical behavior of the corresponding co-crystal as predicted by firstprinciples calculations (Section 5). The computational details are provided in the Appendix.

Methodology
2.1. Density-functional theory for molecules and crystals: One tool for two communities Density-functional theory (DFT) is one of the most successful methods for electronicstructure calculations. Its theoretical foundation lies upon the Hohenberg-Kohn theorems [71], defining the unique correlation between the potential and the density of a many-electron system, which can be used to represent all its ground-state properties. The most popular and efficient recipe to implement DFT is given by the Kohn-Sham (KS) equations [72], which map the many-body problem into an auxiliary system of non-interacting particles ruled by the Schrödinger equation The eigenfunctions of Eq. (1) are the states of the fictitious KS system of independent electrons that are used to express the electron density, and the eigenvalues represent their corresponding single-particle energies. Beside the kinetic-energy term, the Hamiltonian of Eq. (1) includes an effective potential which is given by the sum of three terms: The external potential, v ext , describes the electron-nuclear Coulomb attraction, the Hartree potential, v H [n](r), the mean-field electron-electron coupling, and the exchangecorrelation (xc) potential, v xc [n](r), includes the remaining interactions. While the form of the first two potentials is known, unfortunately, the same is not true for v xc [n](r). Hence, the accuracy of the results provided by the KS equations crucially depends on the approximation chosen for this quantity. The local-density approximation (LDA), originally proposed by Kohn and Sham in their seminal paper from 1965 [72], is based on the homogeneous electron gas model and represents the lowest rank of the so-called Jacob's ladder of densityfunctional approximations [73]. The generalized gradient approximation (GGA) [74] is the next rung, with v xc depending on both the density and its gradients. Higher-level approximations include, on the one hand, the so-called metaGGA functionals [75,76], and, on the other hand, hybrid functionals, where v xc is enhanced by a fraction of Fock exchange [77]. Among the latter, global hybrids are distinguished from the rangeseparated ones, in which the long-and short-range parts of the Coulomb potential are augmented separately [78,79]. Interestingly, the diffusion of these families of functionals has followed separate paths in chemistry and physics, reflecting the substantially different performance of the various approximations on (isolated) molecules and solids. While for inorganic crystals the LDA and GGA have been extensively adopted for decades and specific flavors of range-separated hybrid functionals such as HSE [80] have gained popularity only recently, in the study of molecules the use of hybrid functionals has become a standard since the turn of the turn of the century.
The definition of the basis set for the KS equations (Eq. 1) is another crucial aspect in the implementation of DFT applied to non-periodic and extended systems.
The delocalized character of the electronic wave-functions in solids is prone to a representation in terms of pure plane-waves or to plane-waves augmented by local orbitals [81]. The solution of the KS problem in k-space typically exploits the crystal periodicity and enables the application of efficient parallelization schemes [82]. On the other hand, the non-periodic nature of the electronic states in molecules calls for a realspace description with localized basis sets. Gaussian functions, atom-centered orbitals, and grid representations are the most common options ensuring an optimal trade-off between accuracy and computational costs [83,84,85].

Many-body perturbation theory: Gateway to the excited states
The efficient ab initio scheme provided by DFT is limited by construction to the description of ground-state properties. In order to access excited states, it is necessary to take other routes. One option is to adopt the time-dependent extension of DFT (TDDFT) [86]. This approach is routinely applied both in a perturbative flavor [87] as well as by propagating the KS equations in real time [88]: the former approach gives access to the single-particle composition of the excitations but is limited to the linear response; the latter grants access to the response of all orders [89,90] even in presence of an external, time-dependent electric fields [91,92], but does not provide any information about the nature of the excited states. In both flavors, TDDFT results are crucially dependent of the choice of the exchange-correlation functional. Calculations based on hybrid functionals are largely applied in the chemistry community to describe optical excitations of molecules [69], even in conjunction with embedding methods [93] to account for solvation effects. On the other hand, in spite of enormous efforts in the past decades to develop appropriate approximations [94], TDDFT is not considered reliable to describe excitations in solids, not even when composed of molecules [95,96]. For this reason, its popularity in the solid-state physics community is relatively low.
The state-of-the-art approach to describe electronic and optical excitations in extended systems is many-body perturbation theory (MBPT) [97], including the GW approximation for the self-energy [98] and the solution of the Bethe-Salpeter equation [99]. Modern implementations of MBPT [100,101] are built on top of DFT in the usual frameworks outlined in Section 2.1 for both periodic (see, e.g., Refs. [102,103]) and non-periodic systems [104,105,106]. In this framework, the single-particle Green's function is defined in terms of the non-interacting KS states and energies obtained from Eq. (1). The expression of the dynamically screened Coulomb interaction in addition to the bare Coulomb potential v C , includes the frequency-dependent inverse dielectric function ε −1 , which is evaluated from KS states, too [97]. Eqs. (4) and (5) are the ingredients for the construction of the self-energy in the GW approximation, Σ = iGW . In the original formulation by Lars Hedin [98], this equation was supposed to be solved self-consistently as part of a closed loop of equations. In practice, the perturbative "single-shot" approach G 0 W 0 is typically adopted with good results especially for conventional semiconductors [107]. In this framework, the quasi-particle (QP) equation dressing the KS energies with the self-energy contribution reads [108]: where Z i is a renormalization factor compensating for the evaluation of Σ at KS i , rather than the correct QP i , by means of a linear extrapolation. Starting from the QP-corrected electronic structure, optical excitations are computed by solving the Bethe-Salpeter equation (BSE) [99], the equation of motion for the two-particle correlation function [109]. By construction, this formalism enables the description of excitons quantitatively accounting for electron-hole interactions. In practice, the problem is expressed in terms of an effective two-particle Schrödinger equation where o and u stand for occupied and unoccupied states, respectively. In spinunpolarized systems, the BSE Hamiltonian is expressed as: where the diagonal term,Ĥ diag , accounts for the energy differences between occupied and unoccupied states,Ĥ dir corresponds to the direct Coulomb attraction between the positively-charged hole and the negatively-charge electron, and the exchange termĤ x includes the exchange electron-hole repulsion, Eq. (9) includes the screened Coulomb interaction (Eq. 5), with ε evaluated at ω = 0 (static screening). In Eq. (10),v C is the short-range part of the bare Coulomb potential accounting for local-field effects which are known to be particularly relevant in the optical response of organic materials [95,96,110,111]. The eigenvalues of Eq. (7), E λ , represent excitation energies and mark the resonances in the absorption spectra. In the context of inorganic semiconductors, where excitons are loosely bound and largely delocalized, exciton binding energies, E b , are defined as the difference between excitation energies and the QP optical gap: [112]. In the spectra of these materials, excitons are identified by narrow resonances below the onset [113]. This definition turns out to be inappropriate for organic crystals where the optical absorption spectrum is usually formed by distinct peaks even beyond the absorption onset, which are present already in the (joint) density of states [96,114,115]. In this scenario, a more consistent definition of the exciton binding energy is given by The eigenvectors of Eq. (7), A λ , contain information about the electronic composition and the spatial distribution of the excited states. They act as weighting factors in the transition dipole coefficients, defined as t λ = ou A λ ou o| d|u for localized systems and in terms of the momentum operator, in periodic structures. The imaginary part of the macroscopic dielectric function describes the optical absorption of the material within the unit cell volume Ω. Given the vectorial character of t λ , ε M is a tensor with as many non-vanishing components as enabled by the crystal symmetries. Furthermore, the eigenvectors of Eq. (7) are included in the definition of hole and electron densities, defined as and respectively, for transitions between occupied states φ o and unoccupied states φ u . These quantities were originally introduced for isolated systems [116,117]. In solids, Eqs. (13) and (14) retain the periodicity of the lattice [115]: as such, they highlight the fictitiously periodic distribution of the hole and the electron but cannot provide any information about the exciton localization. For this purpose, it is necessary to introduce the exciton wave function, a six-dimensional quantity expressed by and typically visualized by the isosurface of its square modulus with the hole or the electron position fixed.

Polarizable Continuum Model
The description of electronic structure provided by DFT (and even MBPT) assumes the systems to be in vacuo. Especially for molecules, this is hardly a realistic scenario. Molecules are usually emedded in an environment which has a strong influence on their electronic and optical properties. In order to capture these effects, implicit methods have been developed in quantum chemistry. The polarizable continuum model (PCM) is one of the most popular ones and describes the surrounding medium by a single polarizable continuum interacting with the charge density of the molecule hosted in a cavity. The validity of this approximation depends on the system. Since it fails to capture intermolecular interactions between static dipoles, van der Waals forces, or chemical bonding, including charge transfer, its accuracy is limited to those scenarios in which these effects are negligible, unless corresponding (usually empirical) corrections are included [118].
where v c (r) = 1/|r| is the bare Coulomb interaction and C is the cavity surface, i.e., the interface between the molecule and the medium. In turn, σ[n] is determined from the electrostatic potential created by the molecule's charges on C: where Z j and R j are nuclear charges and positions, respectively. The PCM matrix, Q, encodes all information about the response of dielectric environment. Within the flexible integral equation formalism [119], it is most generally given as where I is the identity operator, and S i , D i and D * i , as well as S e and D e are integral operators, defined by their action on an function u supported on C as [120]: wheren(s) is the cavity surface normal vector and ε is the dielectric constant of the material directly on the exterior of C. The term W (r, r ) appearing in S e and D e is the screened Coulomb interaction in the external medium, disregarding the presence of the cavity. If the so-described environment is homogeneous with dielectric constant ε (a reasonable approximation for solvents or crystalline environments), it is simply given as . In this case, Eq. (18) can be simplified to see Ref. [121]. When instead the medium is anisotropic, adjustments to the formalism are needed, as discussed in Ref. [122] and in Section 4 below.

Organic semiconductors: From molecules to crystals
We start our analysis with the example of anthracene (ANT), the three-ring member of the oligoacene series. In the following, we examine the electronic and optical properties of this system in gas and crystalline phase calculated from DFT and MBPT, and discuss the impact of the chosen modelling framework in reproducing the physics of the system.

Electronic properties
In the first step, we evaluate the frontier states of the isolated ANT molecule from DFT adopting a range-separated hybrid functional and applying the GW approximation in the perturbative approach G 0 W 0 (see computational details in the Appendix). In Fig. 1a), the isosurfaces of the highest occupied and lowest unoccupied molecular orbitals (HOMO and LUMO, respectively) are visualized in a schematic energy-level diagram including the neighboring levels, HOMO-1 and LUMO+1. The picture provided by these results, yielding a fundamental gap of 6.94 eV, is in excellent agreement with the experimental references for the electronic structure of this molecule (6.9 eV, see Ref. [123]). In the bulk phase, ANT crystallizes in a monoclinic structure, space group P 2 1 /a [124], including two inequivalent molecules in the unit cell, see Fig. 1c), where both the real-space representation of the material and its Brillouin zone (BZ) are shown. The resulting band structure, reported in the bottom panel of Fig. 1c), exhibits the hallmarks of organic crystals. First and foremost, the electronic states are no longer described by isolated energy levels but by bands with a finite k-dispersion providing information about the carrier mobility along a certain direction. The relatively flat character of the bands of ANT compared to their counterparts in inorganic semiconductors is related to the molecular nature of the constituents. Since in the energy window displayed in Fig. 1c) the electronic states have π/π * character, the dispersion is larger along those directions where the delocalization of the corresponding electronic cloud is maximized. For this reason, in the electronic structure of longer members of the oligoacene family, such as tetracene and pentacene, the band dispersion is overall enhanced [10,13,52,125]. Another typical characteristic of organic crystals in the band structure of ANT is the twofold splitting of each electronic level, due to the presence of two inequivalent molecules in the unit cell. Depending on the specific direction in reciprocal space, the "sub-states" are either energetically slit up or almost degenerate: the former scenario is realized when the corresponding wave function is spread over both molecules in the unit cell, thereby inducing electronic repulsion; in the latter case, the wave function is localized on only one inequivalent molecule in the unit cell, and thus the mutual repulsion is minimized [52]. As we will see in Subsection 3.2, this characteristic impacts the optical spectrum of the material giving rise to Davydov splitting [52,126]. The computed bandgap, 4.30 eV, indicates a substantial renormalization with respect to result obtained for the isolated molecule (6.94 eV, see Fig. 1a) as a consequence of the screening exerted by the molecules forming the crystal [127]. Notice that the aforementioned value, 4.30 eV, overestimates the experimental band gap of the ANT crystal, ranging between 3.90 and 4.00 eV [128,129,130]. We can ascribe this discrepancy to the fact that the band structure reported in Fig. 1c) is obtained adopting the GGA in DFT, with the QP correction to the band gap included by means of an empirical scissors operator, see details in the Appendix. This approach is often used when experimental references are available (see, e.g., Refs. [11,52,131]), given the high numerical costs of GW calculations on these systems.
In order to understand whether the discussed electronic properties of the ANT crystal are primarily induced by short-or long-range effects, we investigate a nonperiodic cluster formed by the two inequivalent molecules extracted from the unit cell of the crystal (see Fig. 1b). Neglecting the lattice periodicity, the energy levels are obviously dispersionless, confirming the long-range nature of this property. However, the twofold splitting in the energy levels of the molecule is preserved and its magnitude is consistent with the result obtained in the crystal at the high-symmetry point Z (see Fig. 1c). This finding suggests that the short-range interactions between the molecules included in the unit cell are responsible for the Davydov splitting, hinting that this effect can be reproduced even by non-periodic models. The isosurface plots of the orbitals shown in Fig. 1b) reveal the correspondence between the HOMO and the HOMO-1 (LUMO and LUMO+1) of the cluster and the HOMO (LUMO) of the isolated molecule. In both manifolds, the lower-energy state is localized on the same molecule of the cluster. The orbital segregation is consistent with the character of the corresponding states in the crystal at the Γ point. The value of the fundamental gap in the cluster, 6.66 eV, is less than 0.3 eV smaller than the one computed for the isolated molecule. This result, in line with previous findings obtained at the same level of theory for longer oligoacenes [132], is unsurprising and confirms that the band-gap reduction seen in the crystal cannot be captured with such a minimal model. Yet, as discussed in Section 4, a convenient shortcut is available for an accurate evaluation of this quantity at affordable computational costs, without the need to simulating the full periodic crystal.

Optical properties
We now turn to the optical properties of ANT, considering again the isolated molecule, the crystal, and the bimolecular cluster extracted from its unit cell. The absorption spectrum calculated for gas-phase compound (Fig. 2a) is characterized by an optically active excitation at the onset. The energy of this resonance is in excellent agreement with the experimental value [133] and the corresponding excitation is given by the HOMO→LUMO transition polarized along the short molecular axis. This characteristic, common to all oligoacenes [12], can be understood from the parity of the frontier orbitals shown in Fig. 1a). Consistently, the hole and electron densities (Eq. 13 and 14) exhibit the same spatial distribution as the HOMO and the LUMO, respectively. The orientation of the transition dipole moment of this excitation explains its lower oscillator strength with respect to the higher-energy one, polarized along the long molecular axis, which is visible in the spectrum above 5 eV (see Fig. 1a). Moving now to the result obtained for the bimolecular cluster (Fig. 2b), we notice essentially the same spectral features discussed above for the isolated moiety, except that in this case, the lowest-energy peak is formed by two almost degenerate excitations. Their energy separation, of the order of 10 meV in this system, corresponds to the Davydov splitting. The electron and hole densities of the two excitations giving rise to the first peak [see vertical bars in Fig. 2b), left panel] are localized on the same molecule. The lowest-energy excitation stems from the transition between the HOMO and the LUMO+1, while the second one from HOMO-1→LUMO. Their polarization is along the short molecular axis, as in the isolated molecule. Excitations polarized along the long molecular axis appear at higher energies, above 5 eV.
The absorption spectrum calculated for the crystal is plotted showing two diagonal components of the dielectric tensor (see Fig. 2c). A sharp resonance along b dominates the onset followed by a very weak excitation along c, in agreement with earlier results obtained at the same level of theory [11]. The presence of two excitations with two polarizations is due to Davydov splitting. However, in the crystal, this effect is not only related to the presence of two inequivalent molecules in the simulation cell, as in the cluster model, but also to their mutual orientation with respect to the crystalline axes. Appropriate rotation of the dielectric tensor along the electric field directions enhances the spectral intensity of the two Davydov components [52,134]. The spatial distribution of the electron-hole pair corresponding to the lowest-energy exciton of the ANT crystal is shown as the square modulus of exciton wave function, see Eq. (15). In Fig. 2c), the isosurface represents the correlated electron distribution with respect to the fixed hole position, marked in the plot by a red dot. In analogy with its longer oligoacene siblings [12,52], lowest-energy exciton of the ANT crystal has Frenkel character. Although this result may suggest that the adopting a cluster model is sufficient to capture the exciton nature in this system [compare Fig. 2b) and c)], the localization of the electron-hole pair in the crystal is actually a consequence of the large the electron-hole coupling strength, as extensively discussed in previous work on oligoacenes [12,52,60]. Neglecting the proper treatment of these effects, as provided by the solution of the BSE, would lead to an incorrect description not only of the optical spectra but especially of the exciton distribution, which would erroneously retain the periodicity of the electronic wave functions.

Accounting for electrostatic interactions from an implicit environment
In the previous section, we have shown the importance of describing an organic material in its actual state of matter in order to capture correctly its physical properties. Longrange effects such as band dispersion and optical response require accounting explicitly for the lattice periodicity. Modelling the system only through the two (isolated) molecules included in the unit cell renders solely those properties that are ruled by shortrange (i.e., nearest-neighbor) interactions. Yet, the full-fledged simulation of an organic crystal is not a trivial task. Common bottlenecks are the availability of input structures and of the computational resources that enable performing converged calculations on such systems with state-of-the-art ab initio methods. In the following, we see how implicit schemes, accounting only for electrostatic interactions between molecules and their isotropic as well as anisotropic environment and applied in conjunction with DFT, can yield a correct estimate of the band gap of organic materials in different configurations. . This corresponds to the so-called ∆SCF approach, which we apply herein adopting the GGA for the exchange-correlation functional (see computational details in the Appendix).
We start from the gas-phase molecule, for which we obtain IP = 7.07 eV and EA= 0.70 eV, both in good agreement with corresponding experimental data: measurements reported for IP range from 7.40 to 7.45 eV [135,136], while for EA, reference values are between 0.53 and 0.66 eV [137,138,139]. The fundamental gap of gas-phase ANT is therefore equal to 6.37 eV (see Fig. 3a, right), i.e., slightly smaller than the experimental prediction (6.74 -6.92 eV) extracted from the values of EA and IP reported above and reported in Ref. [123]. Discrepancies between calculated and measured values can be ascribed to the adopted GGA for the exchange-correlation functional.
Next, we apply the same procedure embedding the ANT molecule, its anion, and its cation into a PCM cavity filled by a homogeneous medium with ε = 3.10, corresponding to the static dielectric constant computed for the ANT crystal in the random-phase approximation (RPA). In this setup, the calculated IP decreases to 5.79 eV and the EA increases to 1.84 eV (see Fig. 3a, right). It is interesting to notice that both quantities Table 1. Dielectric constants (ε) computed from the RPA for the corresponding crystals, energy gaps (E gap ) of naphthalene (NAP), anthracene, (ANT), and tetracene, (TET) calculated with the ∆SCF method applied to the isolated moieties in the neutral and charged (±1) states supplemented by PCM to mimic the crystalline environment, and experimental reference values: a Ref. [128]; b Ref. [129]; c Ref. [130].  [128,129,130]. It is worth noting that this good match is enabled by the low dispersion of the bands in the ANT crystal (see Fig. 1c), which is in turn an consequence of the localized character of electronic states of the constituents.

System
In the presence of highly dispersive bands, where the energy levels assume very different values at varying k-vector, such a good agreement with an effective model for the environment is not foreseeable. Band dispersion is the consequence of precisely the type of intermolecular interactions that are neglected in the PCM and that require explicit quantum-mechanical modelling to be captured. To prove that the results obtained for ANT are not just a favorable coincidence, we apply the same procedure to other short members of the oligoacene family. The results reported in Table 1 and the shown agreement with experimental references confirm the ability of this approach to correctly reproduce the band gaps of the crystals. It is worth noting that the static dielectric functions used to feed the PCM in these calculations have been computed from the RPA for the corresponding crystalline phases. However, this is not the only way to retrieve these quantities. A cheap alternative is to evaluate them from linear-response calculations on the gas-phase moiety and from the Clausius-Mossotti equation, which connects the microscopic polarizability of the single unit with the macroscopic dielectric constant of the crystal [140,141]. Otherwise, for the considered systems, which are well-known and extensively characterized, corresponding empirical values are available in the literature [129,142].
The generally good accuracy of the adopted ∆SCF approach in the prediction of IP and EA makes viable the usage of optimally tuned, Koopman-compliant hybrid functionals [143]. Parameters therein are adjusted in a non-empirical way until the HOMO energy matches the ∆SCF-determined IP; the EA can be accounted for as well, demanding that it agree with either the LUMO energy of the neutral molecule or -more strictly compliant with Janak's theorem [144] -the HOMO energy of the anion. Actually, the IP and EA obtained with ∆SCF are relatively robust with respect to the replacement of DFT with Fock exchange, meaning that hybrid functionals do not necessarily improve very much the estimate of these values with respect to pure DFT. However, hybrid functionals -particularly the long-range-corrected ones with a large amount of Fock exchange in the asymptotic limit -are essential when computing optical excitations with TDDFT, enabling acceptably accurate results even for chargetransfer excitation energies, which pure DFT notoriously struggles with [145]. Such hybrid TDDFT approaches coupled with the PCM can yield reasonable estimates for bandgaps and even for exciton binding energies [140,141], in spite of neglecting static intermolecular interactions, and artificially confining excitons, e.g., into a single unit cell. These approximations can be softened by replacing a single molecule by a cluster, thus describing a larger portion of the system atomistically and quantum-mechanically. In this case, some pitfalls with respect to symmetry have to be avoided, which is generally lower for clusters than for both isolated molecules and the organic crystals [146].

Repurposing the PCM: From solvent cavities to substrate layers
Using the generic prescription of Eq. (18) for building the PCM matrix, one can model even more complex dielectric environments than an isotropic solvation cavity. The only requirement is an expression for the screened Coulomb interaction W . In Ref. [122], we extended the PCM to describe implicitly an arbitrary stack of layered materials intrinsically characterized by a high degree of anisotropy, introducing LayerPCM. This scheme reproduces the increasingly popular scenario of organic molecules adsorbed on two-dimensional transition-metal dichalcogenides (TMDC) [147,148,149,150,151,152,153]. In this case, the screened interaction in the region above the slab can be expressed as where the so-called polarization contribution, W p , is related to the polarization of the substrate. For a semi-infinite substrate, W p can be written as the electrostatic potential of a single image charge within the substrate region. For a stacked substrate, it is more convenient to calculate it numerically as using a Gauss-Laguerre quadrature [154]. In Eq. (22), ∆r = [(x − x ) 2 + (y − y ) 2 ] 1/2 , J 0 is the zeroth-order Bessel function of first kind [155], and ε(q ) is the effective nonlocal dielectric constant of the substrate, which can be determined using a transfer matrix formalism [122]. An alternative to this Bessel-function-based approach for obtaining W p is the recursive calculation of image charges [156].
We apply LayerPCM to model the hybrid interface formed by an ANT molecule adsorbed on a MoS 2 monolayer (see Fig. 3b). Three parameters are required to model the TMDC (or any other layered material) within LayerPCM: the layer thickness t, as well as the perpendicular and parallel components of the static dielectric constant, ε ⊥ and ε , respectively. We employ a self-consistent scheme to obtain these quantities from first principles. The MoS 2 monolayer is simulated atomistically with periodic boundary conditions in all three directions, inserting a large amount of vacuum in the perpendicular direction to decouple the replicas. The dielectric constants are determined through linear-response calculations within the RPA. It should be noted that corresponding results are effective values, ε ⊥,eff and ε ,eff , describing simultaneously the TMDC and the vacuum layer. In order to extract the dielectric constants of the material excluding the vacuum, we take an initial guess for t and perform the following transformations: where c is the perpendicular lattice constant of the adopted supercell, i.e. the sum of the thickness of the TMDC and the height of the vacuum. The different formulas adopted for the in-plane and out-of-plane directions is due to the different boundary conditions for electric fields with parallel and perpendicular polarization with respect to the dielectric interfaces [157]. Equipped with these values, t is obtained by connecting smoothly the plane-averaged semi-local v xc at short range with the image potential V im of the MoS 2 layer at long range, corresponding to the correct asymptotic tail of the exact exchange-correlation potential [158]. This approach was employed in the past to determine mirror planes for metal surfaces [159]. For the dielectric monolayer, the image potential follows as a special case of Eq. (22) with r = r and ε(q ) corresponding to a single layer: The value obtained for t by smoothly matching v xc and V im is reinserted into Eq. (23) to calculate new dielectric constants. This procedure is iterated self-consistently. For MoS 2 , we obtain ε = 16.64, ε ⊥ = 11.25, and t = 5.4Å. Owing to the prescription in Eq. (23b), ε ⊥ is extremely sensitive to changes in t. Indeed, values for ε ⊥ differing by a factor of 2 have been reported in Ref. [160] using a different definition of t. While such a wide range of results may seem alarming, it can be understood by realizing that the dielectric constant does not single-handedly describe the polarization response of the medium, but the geometry plays a crucial role, too. This is particularly true when the dimensions of the medium are small. Very different combinations of ε ⊥ and t can thus give rise to similar physical scenarios.
Having defined this scheme, we evaluate the band gap of the physisorbed ANT molecule, where the underlying layer is accounted for via its dielectric function only. As in the crystalline environment, the presence of the TMDC reduces significantly the band gap of the moiety. This is a well-known effect that has been reported in a number of first-principles studies on hybrid interfaces formed by conjugated molecules adsorbed on TMDC monolayers [161,162,163,164,165,166]. Comparing the results obtained for isolated ANT and for ANT adsorbed on a MoS 2 monolayer, we notice a reduction of IP and an increase of EA by about 0.8 and 1 eV respectively, which lead to a fundamental gap of 4.57 eV for the physisorbed molecule (see Fig. 3b, right panel). It is worth noting that the difference between IP and EA implies that already the relatively small ANT ions cannot be assimilated to point charges, for which both values would be equal. The band-gap renormalization for the molecule adsorbed on the TMDC is lower than in the crystal. Comparing the results reported in Fig. 3, panels a) and b), this behavior can be ascribed to more significant energetic readjustement of the HOMO level (set equal to -IP, according to Janak's theorem) in the two scenarios. In contrast, the LUMO level, estimated as -EA, differs energetically by 150 meV in the implicit simulation of the crystalline environment and of the monolayer substrate. To analyze more in depth the electronic properties of ANT physisorbed on the MoS 2 monolayer, we perform a fully atomistic benchmark calculation of the whole hybrid system simulating a supercell corresponding to 4×4 the unit cell of the TMDC and with the molecule adsorbed flat at an approximate distance of 3.3Å from the substrate. The electronic structure of the system, evaluated for consistency with the previous runs at the GGA level of DFT (PBE functional, see computational details in the Appendix), features a type-II band alignment between the organic and inorganic components (see Fig. 4). The HOMO of the molecule is situated within the energy gap of MoS 2 , while the LUMO is higher up in the conduction band. We note in passing the clear signatures of mixing between the HOMO-1 of ANT and the wave functions of MoS 2 at the highsymmetry point M in the valence region, which was found to be a hybridization hotspot between TMDCs and polycyclic aromatic hydrocarbons [167]. It is well known that DFT -particularly with semi-local functionals -is not very good at predicting energy level alignments: it notoriously underestimates band gaps as well as molecular IPs and, furthermore, does not capture image-charge induced level renormalization [168,169]. Thus, we correct the frontier energy levels of the hybrid system as well as of the freestanding monolayer for reference with many-body perturbation theory, using the single-shot G 0 W 0 method (Fig. 4b). As expected, this leads to a band-gap opening; the fundamental gap of MoS 2 is increased from 1.70 to 2.76 eV in the hybrid system, and to 2.77 eV in the freestanding configuration, in agreement with previously reported results on a similar level of theory [170,171,172,173]. The presence of the molecule does not lead to a significant band-gap renormalization in the TMDC. The most important change brought about by the QP correction is the inverted order of the occupied energy levels. The severe underestimation of the ionization energy of ANT induced by DFT in the GGA is corrected by a large downshift of the HOMO, while the valence-band maximum of MoS 2 barely changes, resulting in a change from a type-II to type-I alignment with the molecular levels hugging those of the TMDC. Similar observations have been made for other hybrid interfaces [163,165]. Hence, DFT alone yields unreliable results for the alignment of occupied states; for the unoccupied ones, however, the situation is not as dramatic, considering that the GW corrections to the conduction-band minimum and to the LUMO differ only by 0.27 eV, compared to 0.89 eV between the valence-band maximum and the HOMO.
Interestingly, using the ∆SCF+LayerPCM approach described above for calculating the renormalized band gap of ANT, we obtain results in excellent agreement with the GW -corrected ones, even adopting a GGA functional for DFT (Fig. 4b). Taking as a point of comparison the frontier levels of isolated MoS 2 computed from G 0 W 0 , ∆SCF+LayerPCM yields a significantly better estimate for the level alignment of the hybrid interface than DFT alone. The only salient difference between the fully atomistic GW calculation on the hybrid system and the effective treatment is given by a rigid shift of the MoS 2 levels, which is brought about by charge transfer [164]. Since the Fermi energy of ANT is higher than that of MoS 2 , some electronic population is transferred from the organic to the inorganic side of the interface. The resulting partial chargespositive on ANT, negative on MoS 2 -slightly realign the levels of the subsystems with their corresponding electrostatic energy. This effect is obviously not captured when describing the two subsystems individually. However, the amount of charge transferred is based on the level alignment obtained from DFT in the GGA, as we apply only the perturbative G 0 W 0 approach and not the self-consistent flavor of GW . The staggered line-up predicted by GGA can be assumed to lead to an overestimation of the difference between the individual Fermi energies and, thus, of charge transfer. Ongoing work to perturbatively include charge-transfer effects based on the energy levels predicted by calculations on the individual the subsystems is expected to further improve the accuracy of this method.

Donor/acceptor co-crystals
In the last part of this review, we focus on another type of organic materials, namely donor/acceptor complexes and co-crystals. In contrast with monomolecular compounds like ANT, in these systems, a complex interplay between short-and longrange interactions rules the electronic and optical responses [115,174]. In this analysis, we consider the prototypical system formed by quaterthiophene (4T) doped by the electron acceptor tetracianoquinodimethane (TCNQ).
The synthesis and characterization of doped organic films has been the subject of intensive research in the last decade [44,175,176,177,178,179,180]. Depending on the chemical nature of the constituents and their mutual interactions, two main doping mechanisms have been identified, namely integer and partial charge transfer [45]. In the former case, an entire charge carrier is transferred from the donor to the acceptor, leading to efficient (opto)electronic performances. In the latter scenario, the frontier orbital hybridization gives rise to a strong coupling between the components enabling a fractional charge transfer between them. 4T:TCNQ belongs to the second class of systems, also known as charge-transfer complexes (CTC). The synthesis of CTC as ordered co-crystals is still in its infancy. One of the first, successful attempts [181] enabled the resolution of the crystal structure of 4T doped by TCNQ with 1:1 ratio, which is considered in the following.

Electronic properties: Local interfaces vs. long-range effects
4T:TCNQ crystallizes in a triclinic lattice with unit cell containing a donor/acceptor pair with basal planes facing each other [181] (see Fig. 5a). The complex arrangement of the molecules with respect to the crystal structure and their relation with the BZ is schematically visualized in Fig. 5b), where the directions connecting Γ with the other high-symmetry points are highlighted to ease the reading of the band-structure plot (Fig. 5c). By inspecting this result, we immediately notice highest-occupied and lowestunoccupied states sticking out with respect to the manifold of valence and conduction bands, respectively. Interestingly, the character of these states varies at different kpoints in the BZ. At the zone edges (high-symmetry points Z, N , M , and R), the wave functions are spatially segregated over the donor (valence states) and acceptor molecules (conduction states), respectively, see Fig. 5d). Notably, both valence-band maximum and conduction-band minimum are located at M . Elsewhere in the BZ, the electronic states are hybridized between the donor and the acceptor, in analogy with the HOMO and the LUMO of the corresponding molecular complexes [41,44]. This counter-intuitive behavior of the frontier states in the co-crystal can be explained by the long-range nature of its electronic wave functions [115]. To this end, we introduce a tight-binding (TB) model assuming as a basis the molecular orbitals of the isolated bimolecular donor/acceptor cluster extracted from the crystalline unit cell. We assimilate the system to a quantum-mechanical two-level model including the highest-occupied and lowest-unoccupied energy levels, with corresponding wave functions expressed as where C λ m (k) are the orthonormal expansion coefficients for either level λ, and In Eq. (26), ϕ m (r) are the HOMO (m = H) and LUMO (m = L) of the isolated donor/acceptor complex, and R j are the lattice vectors. From Eqs. (25) and (26), it is evident that ψ λk (r) has the same periodicity of the co-crystal and thus fulfills Bloch's theorem. It is worth highlighting the k-point dependence of the expansion coefficients which brings about the k-dependent character of the electronic states within each band seen in Fig. 5c). The wave functions and energies associated to the valence and conduction bands are obtained by diagonalizing the following problem: where the nearest-neighbor (NN) approximation is assumed and the overlap integral is j e −ik·R j ϕ n (r)ϕ m (r + R j )d 3 r ≈ δ mn . Both assumptions are justified by the overall localized nature of the electronic states in molecular crystals. The matrix elements on the left-hand-side of Eq. (27) are calculated as L} are the on-site (R j = 0) and hopping (R j = R a , R b , R c ) integrals. The eigenvectors of Eq. (25) for the wave functions of the valence (VB) and conduction bands (CB) are expressed as and respectively, where with and w(k) = 1 Plugging Eqs. (29) and (30) into Eq.(25), we get: and The energies associated to these states are calculated from Eq. (27) as With this formula, we can finally plot the energies of the VB and CB fit from the TB model (Eq. 36) with those output by the DFT calculations on the co-crystal, see Fig. 6, where we adopt a different k-path with respect to Fig. 5c) to ease visualization. The agreement between the two results is very good and confirms the validity of the model, which, therefore, can be used to rationalize the varying spatial distribution of the VB and CB wave functions in the BZ. To this end, we introduce the so-called segregation factor, and we plot it in Fig. 6 for a direct visualization of the band character. Large values of S(k), found in the vicinity of the high-symmetry points M , R as well as Z and N , correspond to spatial segregation of the wave function on either molecule. On the other hand, the small magnitude of S(k) close to X, Γ, Y , and L is indicative of orbital delocalization across the donor/acceptor interface. We can rationalize these results going back to the equations introduced above. From Eq. (31), it appears that µ(k) is maximized when A HL (k) is large, i.e., when there is nonnegligible mixing between the HOMO and the LUMO in the donor/acceptor complex. This happens along the crystallographic orientations that are substantially collinear to the lattice vectors aligned with the stacking direction of the molecules (see Fig. 5b). Since the off-diagonal elements ofĤ T B are non-negligible, the expansion coefficients of the wave-functions in Eqs. (29)- (30) have similar magnitudes. In this scenario, longrange interactions prevail over local couplings between donor and acceptor moieties, and induce the spatial segregation of the frontier states seen in Fig. 5d). On the other hand, along the directions that are approximately parallel to the molecular planes and thus orthogonal to the stacking direction, A HL (k) → 0 and no mixing between HOMO and LUMO takes place. In this case,Ĥ T B becomes diagonal, and one C λ m (k) coefficient dominates over the others. Consequently, µ(k) → 0 implying S(k) → 0: local couplings prevail over the long-range interactions, and the hybridized character of the frontier states is preserved in the periodic wave functions (Fig. 5d).

Optical Properties: Charge-transfer excitations enhanced
With the gained understanding of the electronic properties of the 4T:TCNQ co-crystal, we now analyze its optical excitations. Similar to ANT (Section 3 and Ref. [11]) and to organic crystals in general [7,10,52,96,110,111,13], the absorption spectrum of 4T:TCNQ is dominated by a strong resonance at the onset (see Fig. 7), on which we focus in the following. The pronounced anisotropy of the co-crystal is apparent in both the real and imaginary part of the dielectric tensor computed from the BSE [Fig. 7a), left and middle panels]. All the six inequivalent components contribute to the lowest-energy excitation. The one with the largest spectral strength is along zz, which is roughly aligned to the stacking direction of the molecules in the unit cell [see Fig. 5a)-b)]. This behavior is consistent with the first optical excitation of the isolated donor/acceptor complex [115], visualized in the spectrum in Fig. 7a), right panel. The main contribution to the lowest-energy excitation in the co-crystal comes from the transition between the VB maximum and the CB minimum, again in line with the behavior of the cluster where the first excited state corresponds to the HOMO→LUMO transition [115]. However, as discussed above, the wave functions of the frontier states are segregated on the donor and acceptor molecules, respectively (see Fig. 5d). As a result, the first exciton in the co-crystal has charge-transfer character, as indicated in Fig. 7b) by the plots of the hole and electron densities, opposite to the Frenkel nature of its counterpart in the isolated complex (see Fig. 7c). In light of this behavior, the strong intensity of first excitation of the co-crystal may appear puzzling. To solve this conundrum, we make use again of the TB model introduced in Section 5.1.
Starting from the definition of Bloch's states (Eq. 26), we peel off the periodic part and use it to express the matrix elements for the optical transitions as is the transition dipole moment between the (localized) molecular orbitals of the donor/acceptor cluster extracted from the unit cell of the co-crystal. This term is not the only one entering the expression of the matrix elements. The first term on the right-hand side of Eq. (39), including the k-derivative of the expansion coefficients introduced in Eqs. (29)-(30), enables long-range couplings between the periodic states even when their spatial overlap is small. Thanks to this contribution, charge-transfer excitations like the lowest-energy one in the spectrum of the 4T:TCNQ co-crystal and of its siblings with partially or fully fluorinated acceptors [115] exhibit finite oscillator strength. Notice that, in contrast, charge-transfer excitations in non-periodic donor/acceptor organic interfaces are characterized by very weak spectral intensity [179,180], due to the negligible spatial overlap of the wave functions involved in the transition.

Summary and conclusions
We have presented an extended overview on the electronic and optical properties of different kinds of organic materials computed from first principles. Adopting a solidstate physicists' perspective, we have described the systems in various states of matter, ranging from isolated molecules in gas-phase and in solution to extended crystals. With the example of anthracene, we have shown that accounting for the native periodicity of the material is particularly relevant in order to capture the k-dispersion of the electronic bands and to provide a reliable starting point for the calculation of the optical properties, including excitons. On the other hand, mimicking the crystalline environment by immersing a single molecule in a uniform electrostatic cavity, as enabled by the PCM, leads to very accurate values for the fundamental gap as well as for electron affinity and ionization potential. This result suggests that when these properties are sought, an effective approach, considerably cheaper in terms of computational efforts than the full-fledged simulations of the whole crystal, can be successfully employed. Similarly, adopting LayerPCM, the extension of PCM to deal with anisotropic layered substrates that we recently introduced [122], grants access to the correct band-gap renormalization of conjugated molecules adsorbed on a TMDC monolayer. The application of this method to other low-dimensional, polarizable substrates is straightforward. Finally, as the last example, we considered a a charge-transfer complex formed by p-doped oligothiophene and inspected the mutual influence of short-and long-range interactions on its electronic and optical properties by simulating the system both as an isolated cluster and as a periodic co-crystal. We found that only accounting for the lattice periodicity, one can correctly reproduce the intrinsic Bloch character of the wave function in the lattice, which crucially impact on the electronic and optical properties of the system. These effects can be rationalized with the aid of a tight-binding model in the basis of the donor/acceptor unit. This approach not only unfolds the physics of the problem in a transparent manner but, most importantly, draws a connection to consolidated theoretical methods for the description, for example, of the transport properties in this class of systems [182,183,184]. In our analysis, we have focused on the strengths of established solid-state physics approaches such as DFT, alone and coupled with effective models such as PCM and LayerPCM, and MBPT. Of course, there are many aspects related to the properties of organic materials that cannot be quantitatively captured with these methods. For example, exciton dissociation and the subsequent dynamics of the photogenerated charge-carriers can be hardly described at the mean-field level of DFT and even MBPT cannot directly deliver more information than exciton distribution and binding energies. Model Hamiltonians accounting for all relevant interactions are much for suitable to accomplish these tasks [59,185,186]. Likewise, intersections with machine learning have led to steps forward in crystal structure prediction of organic materials [187] even with targeted properties, such as enhanced singlet fission [188]. These few examples shows the importance of establishing a common ground for truly interdisciplinary research in the field of organic materials. By offering our view on this topic, we hope to have provided our contribution in this direction.

Data Availability Statement
The data that support the findings of this study are openly available on Zenodo at the following links: 10.5281/zenodo.7334160 (electronic and optical properties of anthracene molecules crystal); 10.5281/zenodo.7333310 (electronic and optical properties of donor/acceptor co-crystals); 10.5281/zenodo.7329775 (molecules on implicit substrates).

Anthracene molecule and crystal
The results reported in Section 3 for anthracene molecules and crystals are obtained through a two-step procedure. First, all structures are optimized by minimization of the interatomic forces until the threshold of 10 −3 eV/Å. In the optimization of the molecular clusters, the positions of C atoms are clamped in order to preserve the herringbone angle between the two molecules assigned in the input structure to mimic their arrangement in the crystal. Otherwise, the dimer relaxes to a configuration with the basal planes of the two molecules facing each other. These calculations are performed with the all-electron code FHI-aims [84] adopting tight integration grids, TIER2 basis sets [189], and the generalized-gradient approximation in the Perdew-Burke-Ernzerhof (PBE) parameterization [190] for the exchange-correlation potential. Van der Waals interactions are included via the pairwise Tkatchenko-Scheffler scheme [191].
In the second step, electronic and optical properties are computed for the obtained geometries. For the finite systems (isolated molecules and clusters), the MOLGW code [106] is employed with Gaussian-type cc-pVTZ basis sets [192] in conjunction with the frozen-core approximation and the resolution-of-identity approximation [193]. The range-separated hybrid functional CAM-B3LYP [194] is adopted in the DFT step as a starting point for the subsequent G 0 W 0 and BSE runs. Calculations on the crystal are performed using the all-electron full-potential code exciting, implementing the family of linearized augmented plane-wave plus local orbitals (LAPW+LO) methods [81]. Muffin-tin radii of 1.05 bohr and 0.7 bohr are chosen for carbon and hydrogen, respectively, along with a plane-wave cutoff G max =4.7 bohr −1 . The GGA in the PBE parameterization [190] is used for the exchange-correlation potential. In the DFT calculations, the BZ is sampled by a 8 × 6 × 5 k-mesh, while in the BSE step [103] a Γ-shifted 8 × 6 × 5 k-point mesh is adopted. The QP correction is added to the KS gap through a scissors operator equal to δ =2.27 eV, extimated by aligning the lowestenergy excitation computed for the ANT crystal with available experimental data [195]. In the band-structure plot shown in Fig. 1c), the binding energy of the first exciton, E b =1.30 eV, output by our BSE calculation, has been added to δ, in analogy with the procedure adopted in Ref. [52]. In the solution of the BSE, the screened Coulomb potential W is computed using 200 empty bands and the energy threshold for the local field effects is set to 1.0 Ha. In the construction and diagonalization of the BSE Hamiltonian, 10 occupied and 10 unoccupied bands are considered.

DFT+PCM and LayerPCM
DFT calculations coupled with the PCM are performed with a locally modified copy of the Octopus code (v9.2) [196]. Nuclear potentials and core electrons are approximated by norm-conserving SG15 pseudopotentials [197], the exchange-correlation potential by PBE [190]. Neutral and singly-ionized molecules are treated in spin-restricted and unrestricted frameworks, respectively. KS wave functions are represented on a realspace grid generated by uniform sampling with spacing 0.21Å of a union of spheres of radius 4.5Å centered at each atomic position. The PCM cavity is generated in a similar, yet not identical fashion [198], using smaller spheres with atom-dependent radii. Here, we take the van der Waals radius of the respective element, scaled by a factor of 1.1.
For the atomistic calculation of the isolated TMDC monolayer and the ANT@MoS 2 interface as well as for the determination of the dielectric constants fed into the PCM, we exploit the interface between the Quantum Espresso suite (v6.8) [199] and the Yambo code (v5.0.2) [102,200]. With the former, we obtain the DFT starting point for many-body calculations performed with the latter. We again employ SG15 pseudopotentials [197] for the cores and PBE for the exchange-correlation functional [190], neglecting spin-orbit coupling. Structural relaxations are done with a wave-function cutoff of 60 Ry and including the pairwise Tkatschenko-Scheffler dispersion correction [191], reducing the components of all forces to below 10 −3 Ha/bohr. At this step, we use k-point samplings of 8 × 8 × 1, 2 × 2 × 1, and 4 × 6 × 4 for MoS 2 , ANT@MoS 2 , and the acene crystals, respectively. To simulate the ANT@MoS 2 interface, we choose a lattice constant of 20Å in the non-periodic direction including a sufficiently large vacuum region to decouple the replicas, and we include a dipole correction. We use Wannier interpolation [201] as implemented in the Wannier90 code (v3.1.0) [202] to determine the density of states of ANT@MoS 2 on a 100 × 100 × 1 k-grid, allowing for the application of only minimal artificial broadening in the electronic structure, which endows the obtained linewidths with physical meaningfulness. The corresponding band structure is unfolded to the BZ of the MoS 2 unit cell [167,203]. The DFT calculations for the generation of the required large number of unoccupied orbitals in the many-body simulations are done with a wavefunction cutoff of 40 Ry, sticking with the same 4 × 6 × 4 k-mesh for the acenes, but increasing it to 24 × 24 × 1 and 6 × 6 × 1 for MoS 2 and ANT@MoS 2 , respectively, which is necessary to obtain converged GW results. We calculate the dielectric functions of the acenes by means of linear-response calculation on the level of the RPA, featuring 300 conduction bands and a G, G sum kinetic cutoff of 3 Ry. The values used in the PCM correspond to the trace of the ω = 0 dielectric tensor, averaging xx, yy, and zz components. For MoS 2 , we take 100 conduction bands and a cutoff of 5 Ry to calculate the (effective) anisotropic dielectric functions within the RPA.
The GW calculations on MoS 2 and ANT@MoS 2 are performed with the Yambo code [102,200] adopting the plasmon-pole approximation for the dynamically screened interaction W . 300 (MoS 2 ) and 500 (ANT@MoS 2 ) conduction bands are included for both W and the correlation part of the self-energy, Σ c and employ the sum-over-states terminator [204]. The G, G sum is cut off at kinetic energies of 5 Ry. A Coulomb cutoff [205] along the out-of-plane direction is applied, truncating the interaction at 19.57Å (c = 20Å). Coulomb integrals over the BZ for the smallest ∼300 G vectors are calculated with the random integration method, replacing the sum over a uniformly sampled Brillouin zone by a 10 6 q-point Monte-Carlo integration. Here, we neglect the anisotropy of W for small q, choosing the corresponding polarization vector as (1,1,1) √ 3 , representing an orientational average.

Donor/acceptor complex and co-crystal
The results for the 4T:TCNQ co-crystal presented in Section 5 have been obtained from DFT using the code Quantum Espresso [199] and from MBPT using Yambo [102,200]. In the former step, a plane-wave basis set cutoff of 50 Ry (200 Ry) for the wavefunctions (electron density) and norm conserving pseudopotentials [206] are employed. The PBE functional [190] is adopted in conjunction with the Tkatchenko-Scheffler pairwise scheme [191] to account for van der Waals interactions. A 4×4×4 k-grid is used to sample the BZ. Self-consistent calculations are carried out on optimized geometries obtained by minimizing interatomic forces with a threshold of 10 −4 Ry/Å. BSE calculations are performed on top of the DFT electronic structure with the quasiparticle correction added in the form of a scissors operator of 1.41 eV to the conduction bands based on available experimental spectra [44]. The BSE Hamiltonian is constructed including 20 occupied and 40 unoccupied states, and by 4×4×4 k-point grid; it is diagonalized using the Haydock-Lanczos algorithm [207].
Calculations on the 4T:TCNQ cluster are performed with MOLGW [106] adopting the same computational parameters employed for ANT, except for the Gaussian basis set, which is here reduced to aug-cc-pVDZ.