Accurate description of charged excitations in molecular solids from embedded many-body perturbation theory

We present a novel hybrid quantum/classical (QM/MM) approach to the calculation of charged excitations in molecular solids based on the many-body Green's function $GW$ formalism. Molecules described at the $GW$ level are embedded into the crystalline environment modeled with an accurate classical polarizable scheme. This allows the calculation of electron addition and removal energies in the bulk and at crystal surfaces where charged excitations are probed in photoelectron experiments. By considering the paradigmatic case of pentacene and perfluoropentacene crystals, we discuss the different contributions from intermolecular interactions to electronic energy levels, distinguishing between polarization, which is accounted for combining quantum and classical polarizabilities, and crystal field effects, that can impact energy levels by up to $\pm0.6$ eV. After introducing band dispersion, we achieve quantitative agreement (within 0.2 eV) on the ionization potential and electron affinity measured at pentacene and perfluoropentacene crystal surfaces characterized by standing molecules.


I. INTRODUCTION
The ability to accurately predict the energies of charged excitation from first principles is of primary importance for the computational study of organic conjugated materials that find applications in electronic and opto-electronic devices, since phenomena such as charge injection from a metal electrode, electron-hole separation in solar cells or their radiative recombination in light-emitting diodes do all crucially depend on the energetics of electronic energy levels. [1][2][3][4] The achievement of quantitative accuracy on quantities such as the ionization potential (IP) or the electron affinity (EA) is definitely a challenging task for theory. Well-grounded approximations and efficient implementations are both required to make calculations on systems counting a large number of atoms accurate and feasible at the same time. An additional hurdle comes from the subtle effect of intermolecular interactions in the solid state. In fact, while IP and EA are well-defined properties of a given (isolated) molecule, the same quantities can present variations that can exceed 0.6 eV between different solid samples of the same compound. [5][6][7] Such a large variability originates from intermolecular interactions of electrostatic nature and reflects in charged excitations energetics that depend on morphology and, in crystals, on the facet though which electron are injected or extracted. 5 IP and EA are therefore not intrinsic properties of a given compound but instead depend on the molecular organization, e.g. amorphous vs. crystal or standing vs. lying molecules. 3,4 Many-body perturbation theory (MBPT) techniques, such as the Green's function GW formalism, 8,9 stand as the state-of-the-art for the first principles description of charged (quasiparticle) excitations in condensed matter. Originally developed within the solid state physics community, the GW formalism has been extensively applied to inorganic solids in the last decades, 10,11 leading to a substantial improvement over density functional theory (DFT) in the description of the electronic band gap. The GW formalism is gaining increasing attention also in the context of organic systems, with several applications to extended solids reported recently using periodic boundary condition implementations, [12][13][14][15][16] including studies of the band structure of the prototypical molecular semiconductor pentacene (PEN) [12][13][14] and subsequent investigations of optical properties within the Bethe-Salpeter formalism [12][13][14][16][17][18][19][20] .
A severe limitation of periodic bulk calculations is that they cannot attain the absolute value of charged excitations due to the missing internal reference for the energy of a freeelectron. The use of slab geometries allows the definition of a consistent vacuum level and hence the calculation of absolute values for IP and EA. However, such calculations come at a high computational cost and are difficult to converge with respect to slab and vacuum thickness, making this route impractical especially for organic solids with many atoms in the unit cell. Kang et al. 21 very recently proposed a strategy to obtain IP and EA combining bulk GW calculations with DFT slab calculations, the latter permitting to refer GW quasiparticle energies to the appropriate surface-specific vacuum level. 21 Such an approach is, however, loosely consistent since it does not account for the reduced screening at crystal surfaces, hence results in a systematic underestimation of the gap. 21 The development of GW implementations on Gaussian atomic orbital bases greatly facilitated GW calculations of finite, aperiodic, systems. [22][23][24][25][26] The GW formalism achieved a very accurate description of quasiparticle energies and gap in isolated molecules, as demonstrated by extensive benchmarks against gas-phase experiments 22,[27][28][29][30][31][32][33] and high-level quantum chemistry calculations. [34][35][36][37][38][39] Thanks to efficient algorithms and parallel implementations, GW calculations enabled accurate calculations on systems exceeding hundred atoms. 23,[40][41][42][43][44] Yet, charged excitations in extended systems are largely governed by long-range electrostatic interactions that are not amenable to a full QM treatment.
In this paper, we present a novel hybrid quantum/classical (QM/MM) approach to quasiparticle excitations combining a state-of-the-art implementation of the GW formalism 22 for the QM subsystem with a discrete polarizable model of atomistic resolution. Such an approach goes beyond pioneering implementations that describe the MM dielectric medium as a regular grid of polarizable centers 45 or with the polarizable continuum model. 46 In the present work the MM subsystem is described by the charge response (CR) model by Tsiper and Soos 47 that provides a careful description of the static dielectric response of molecular solids. [48][49][50] The CR and related microelectrostatic models 51,52 greatly contributed to the comprehension of the role of intermolecular electrostatic and polarization interactions on photoelectron spectroscopy measurements 53,54 and on the energetics of charge carriers in organic solar cells. [55][56][57] Our hybrid formalism has been very recently applied to the pristine 58 and doped 42 PEN crystal. A first validation of our hybrid scheme has been demonstrated by comparing our results for the electronic band gap of PEN to experimental data 6,59 and values provided by plane-wave GW calculations. 13 In the present work the embedded GW framework is extended to account for electrostatic crystal field effects, hence providing access not only to the band gap, but also the absolute values of IP and EA at specific crystal surfaces. Our hybrid methodology is applied to model photoemission spectra of PEN and perfluoropentacene (PFP), two widely studied molecular semiconductors for which accurate photoemission data are available for solids of well-defined surface structure. 6,59 The quantitative (within 0.2 eV) agreement on IP and EA for both compounds in the gas and solid phase demonstrates the accuracy and the internal consistency of our approach.
The paper is organized as follows. Section II describes our hybrid formalism in full detail.
Results for PEN and PFP are presented in Section III, where we highlight the importance of crystal field effects on the charged excitations of individual embedded molecules in the bulk and at crystal surface. Our results are discussed and compared to experiments in Section IV, where we dissect the different contributions from intermolecular interactions to IP and EA in PEN and PFP solids, i.e. polarization, crystal field and band dispersion, the latter being accounted for with an ab initio parametrized tight binding model. The main conclusions are finally drawn in Section V.

II. QM/MM METHODOLOGY
As in other hybrid QM/MM approaches, our embedded MBPT calculations are defined by the level of theory employed for the QM and MM subsystems and by the formalization of the interaction between the two parts. We provide below the details of the theory along with approximations and expedients employed to make our hybrid framework feasible and computationally efficient.
A. MM embedding of the ground-state DFT calculation The anisotropic charge densities of neutral organic molecules are source of intense and inhomogeneous electric fields within the crystals that can affect IP or EA by several tenths of an eV, as summarized in recent review papers. 3,4 We account for crystal field effects at the DFT level (providing the starting point for the subsequent GW treatment) by computing the Kohn-Sham (KS) eigenstates {φ n } and eigenvalues {ε n } of the QM subsystem in the field of the surrounding neutral molecules in the crystal, described at classical MM level.
The MM embedding implies a modification of the electrostatic potential experienced by the QM system, i.e. where the charged excitation is created. Such an electrostatic effect is well described in a ground-state DFT calculation and should not be confused with the dynamic reaction of the system to the ionization, which is accounted for within the GW formalism (see Section II C).
Such a strategy, namely starting MBPT calculations with DFT eigenstates obtained in the electric field of the classical environment, was also recently applied for the study of the optical properties of molecular systems in condensed phases. 60-62 As shown below, the effect on the absolute position of the occupied/virtual electronic energy levels is significantly larger than that on the optical excitations relying on energy difference.
An accurate classical description of the MM subsystem can be obtained with discrete polarizable models of atomistic resolution describing the static (zero-frequency) dielectric response of systems of interacting molecules. Our calculations can in principle account for different contributions (i.e. ionic, vibrational, electronic) to the dielectric susceptibility, yet, in the present case of organic crystals of neutral molecules, only the leading electronic response is considered. 49,50 Specifically, in the present work we resort to the charge response (CR) model by Tsiper and Soos,47 describing the anisotropic linear molecular response to electric fields in terms of induced atomic charges and induced dipoles. The CR model is entirely parametrized with quantum-chemical calculations. It has been shown that CR models, including the very similar charge response kernel theory by Morita and Kato, 63 provide a quantitative description of the static permittivity tensor of several molecular crystals. [48][49][50]64 A careful description of the electrostatic potential of isolated neutral molecules is of particular importance for an accurate assessment of crystal field effects within the MM model. In our CR scheme we rely on point atomic charges obtained from the fitting of the electrostatic potential generated by the DFT electron density. 65 From a practical point of view, an iterative scheme consisting of cross-coupled DFT and CR calculations is set-up to obtain KS orbitals in the self-consistent field of permanent and induced multipoles in the MM region. We start from a gas-phase DFT calculation on the QM subsystem and compute the electric potential and fields generated by the DFT electron density at the atomic sites of the MM region using efficient and accurate Coulomb-fitting resolution-of-the-identity (RI-V) techniques. Fields and potentials generated by the QM electron density are then used to compute the induced charges and dipoles within the MM region, accounting for mutual interactions between MM molecules. The DFT calculation for the QM subsystem is then repeated in the field of permanent and induced multipoles in the MM region, and the whole procedure is iterated until achieving self-consistency. For the crystalline materials considered in this work, the energies of occupied and virtual molecular orbitals converge (within 1 meV) in 3 iterations.
Electrostatic interactions are notably long-ranged and special care is required when truncated sums are employed to approximate results for infinite systems. 4,66 Moreover, the potential generated by ordered arrays of quadrupolar molecules does not only depend on its size but also on its shape, leading to different values in crystalline bulk (3D geometry) and in 2D slabs, the latter depending on the crystallographic facet exposed to the vacuum. Results We describe briefly the GW formalism on which hinges the chosen QM framework within the QM/MM approach developed in this study, mostly emphasizing the main features related to the embedding strategy. More details about MBPT can be found in review articles devoted to the GW approach. [8][9][10][11][67][68][69] Our starting point is the time-ordered one-body Green's function G describing the propagation in time of an added (removed) electron to (from) the N -electron system in its ground-state. More precisely, G reads: where ψ GS (N ) is the N -electron ground-state wave function and ψ (r, t),ψ † (r, t) are the destruction/creation field-operators in the Heisenberg representation. Alternatively, it can be shown that G adopts in a frequency representation the form: To proceed further and obtain the G operator in practice, one should solve the following equation-of-motion: where the one-body Hamiltonian h 0 contains the kinetic energy operator and the ionic and Hartree potential. The self-energy operator Σ(r, r ; ω) represents all exchange and correlation effects. Note that it is non-local and energy dependent, in contrast e.g. with adiabatic and semi-local DFT exchange-correlation functionals. While such formal developments are exact, an accurate approximation for the self-energy Σ is required. Within the GW formalism, which can be considered as the lowest-order approximation to Σ in terms of the screened Coulomb potential W , 8,67 the quantities to be calculated read: where v(r, r ) = (r − r ) −1 is the bare Coulomb potential, χ 0 the independent-electron susceptibility and W the screened Coulomb potential. The {f n } are occupation numbers and {ε n } are the energies of the KS eigenstates (orbitals) {φ n } that will be corrected within the present GW formalism.
While the Green's function G can be calculated with the above set of equations, in practice, and since Σ contains all effects related to exchange and correlation, the most common approach consists in replacing in a perturbative fashion the DFT exchange-correlation potential V DF T XC contribution to the KS eigenstates by its self-energy analog, namely: Starting from the KS eigenstates obtained with a given exchange and correlation functional, Σ GW is constructed following Eqs. 3-6. The GW quasiparticle excitations E GW n are then obtained by correcting the KS energy levels according to Eq. 7. Such a scheme is labeled G 0 W 0 , where the "0" subscript indicates that G and the screened-Coulomb potential W are built from the zero-ordered (uncorrected) KS eigenstates. The G 0 W 0 scheme provides improved quasiparticle excitations with respect to DFT, leading, however, to results that depend on the starting exchange and correlation functional. 35,37,39,71,72 A more accurate, although computationally more expensive, approach consists in achieving a partial self-consistency on the eigenvalues only. In the so-called evGW approach, the many-body corrected energies are in fact re-injected in Eqs. 3-6 (E GW n → ε n ) building the self-energy operator corrected to the next order, and the whole procedure is iterated until convergence of {E GW n }. The dependence on the starting functional is significantly reduced in the evGW scheme, leading to quasiparticle excitations in quantitative agreement with experimental values or higher level CCSD(T) calculations. 35,36,38 C. MM dielectric contribution to the screened Coulomb potential W We now show how the contribution of the MM subsystem dielectric response can be merged with the GW formalism to properly contribute to the energy required for adding or removing an electron to the QM subsystem. On general grounds, the analysis of Eq. 6 shows that if two subsystems, hereafter labeled "1" and "2", have non-overlapping orbitals, the independent-electron polarizability cannot couple the two systems, namely χ 0 (r 1 , r 2 ) is zero for any pair of positions r 1 and r 2 in systems 1 and 2, respectively. As a consequence, the screened Coulomb potential restricted to the QM subsystem (1), W 11 , reads: using a block notation where index 1 (2) corresponds to points located in area 1 (2). After some algebra, one obtains the following set of equations: whereṽ 11 is the Coulomb potential in the QM cavity screened by the MM subsystem only (namely without the response of the QM section itself that is incorporated in the χ 0 11 susceptibility in Eq. 10) and χ * 22 is the interacting polarizability of system 2 alone, i.e. in the absence of system 1. Upon introducing real-space coordinates, the Coulomb potential within the QM region is renormalized by adding v reac (r 1 , representing the reaction field generated in r 1 by the MM subsystem in response to a charge added in r 1 , with both r 1 and r 1 pointing in the QM subsystem 1. The reaction field v reac in Eq. 13 is therefore the key quantity through which classical polarizabilities of molecules belonging to the MM region enter the embedded GW calculation restricted to the QM subsystem. Such a calculation can then be performed as a standard GW calculation in the gas phase, but with the bare Coulomb potential substituted by the renormalized (MM-screened) potentialṽ 11 . The construction of the GW self-energy actually requires the knowledge of the dynamically screened W (ω) Coulomb potential in Eq. 5, accounting for the fact that the system dielectric response is frequency-dependent in the optical range. Indeed, in a full GW calculation the dispersion of the QM subsystem susceptibility is accounted for in Eq. 6, while a possible frequency dependence of the MM subsystem polarizabilities would result in a frequency-dependent reaction field.
Although the dependence of the optical dielectric properties on photon frequency is experimentally well documented in organic solids, 73 the classical polarizable models we rely on focus on reproducing correctly the correct optical dielectric response in the low-frequency, ω → 0, limit. While generalizing MM polarizable models to dynamical response may be considered, we describe here an alternative strategy that consists in the merging of static MM polarizable models with the static limit of the GW formalism. Such an approach has been recently applied to couple the GW formalism to continuum 46 and discrete polarizable models. 58 The static formulation of the GW formalism, the so-called static Coulomb-hole plus screened exchange (COHSEX) approximation, was discussed in the seminal paper by Hedin 8 and was shown to be very efficient within the framework of simplified self-consistent GW calculations. 74 COHSEX calculations are known to be less accurate than the full GW ones in determining e.g. the band gap of semiconductors. However, within the purpose of the present QM/MM scheme, the COHSEX approximation will only be adopted to obtain the contribution of the MM environment to quasiparticle excitations. Using a symbolic notation, we decompose the self-energy operator for the embedded QM system as follows: This formula approximates the self-energy operator of the embedded system by the sum of its analogue for the isolated QM system, plus a correction calculated at the COHSEX level, namely as the difference between the self-energy obtained including or not the MM reaction field contribution. The reason for such a formulation is that the use of the COHSEX approximation in the form of a difference allows to reduce the error introduced by replacing the frequency-dependent optical dielectric constant by its low-frequency limit. Even though the screening potential W in the quantum-mechanical region is modified by the MM response, one here assumes that the dynamical screening contribution entering in the difference (Σ GW/MM − Σ COHSEX/MM ) largely cancels with the one in (Σ GW − Σ COHSEX ).
The approximated GW /MM self-energy in Eq. 14 can be finally used to compute the quasiparticle energies basis is lower than 10 meV in the case of PEN. 58 Both the long-range nature of the reactionfield, and the fact that polarization energy is calculated as an energy difference, may explain this observation.

D. Reaction field matrix on a Gaussian basis
In our Gaussian atomic orbital implementation we do not calculate v reac on a (r, r ) real-space grid but look for the following matrix elements: v reac (β, β ) = dr dr β(r)v reac (r, r )β (r ) (16) namely the two-center two-electron Coulomb integrals between auxiliary Gaussian orbitals {β} located at the atomic sites in the QM region. This auxiliary basis stems from the Coulomb-fitting resolution-of-the-identity (RI-V) formulation of the GW implementation we adopt.
In practice, before performing the GW calculation, we compute the self-consistent rearrangement of MM charges and dipoles induced by the potential generated by the charge density associated to each orbital β of the auxiliary basis, namely where Since the {β} orbitals serve as a basis set to represent charge densities, the reaction field Another strategy is the extrapolation of the reaction field matrix elements v reac (β, β ) to the infinite MM cluster limit before performing a single COHSEX calculation that will directly target the infinite crystal. In the case of auxiliary s orbitals (angular momentum l = 0), representing electrical monopoles in the RI-V formalism, the reaction field matrix scales as R −1 in three dimensions. For β (β ) orbitals of arbitrary angular momentum, the reaction field matrix elements scale as R −(1+l+l ) , allowing straightforward extrapolation. 75 The decay of v reac (β, β ) elements is faster and faster for high l, leading in some cases to values that are practically already converged in clusters of relatively small size. Additional details on the efficient calculation of the reaction field matrix are given in Appendix A.
Once the extrapolated reaction field is obtained, the polarization energy can be directly obtained in the infinite cluster limit by performing a single COHSEX calculation. The result obtained with this approach for the HOMO and LUMO polarization energies of PEN (filled dots in Figure 1) are practically identical to those extrapolated in the first brute force scheme.
In the following we will describe both calculations of charged excitations in a bulk material and at its surface. The reaction field matrix for surface calculations is extrapolated in the limit of an infinite semi-sphere, strictly describing the polarization response of a semi-infinite crystal to the charging of a molecules at the surface. Such an approach can be considered a good approximation also for molecular films on insulating substrates of comparable dielectric constant, such as SiO x .

III. CHARGED EXCITATIONS IN PENTACENE AND PERFLUOROPENTACENE
Our novel Green's function QM/MM formalism is here applied to calculate charged excitations in PEN and PFP crystals, considering the ionization of molecules in the bulk and at the crystal surface. The photoelectron spectra of PEN and PFP have been extensively studied experimentally 6,59,86 and theoretically, 3,54,59 as they represent an ideal case study for dissecting the different contributions to charge transport levels arising from intermolecular interactions in the solid state. 4 Photoelectron spectra of PEN and PFP are therefore chosen here to demonstrate the accuracy and internal consistency of our hybrid formalism.
The two molecules have indeed similar chemical and electronic structure (the elementary Hückel model for π-electrons cannot distinguish between the two), leading to comparable frontier orbitals 87,88 and nearly identical polarizability tensors. 54 Both compounds crystallize in a layered structure characterized by planes of nearly standing (slightly tilted) molecules arranged in a herringbone fashion. Such a layers are parallel to the (001) and (100) planes in PEN and PFP, respectively, which are low energy crystal facets usually found also in molecular films grown on insulating substrates such as SiO x . 6 A crucial difference between the two molecules resides in their ground-state electrical quadrupole moment, which has principal components of comparable magnitude but opposite sign, owing to the different polarity of C-H vs. C-F bonds. 54,89 The availability of accurate experimental data from ultraviolet photoemission spectroscopy (UPS) and low-energy inverse photoemission spectroscopy (LEIPS) for crystalline films of standing molecules 6,59 make these systems ideal for benchmarking purpose. UPS and LEIPS spectra have been also reported for films of laying molecules on metal or graphite substrates. 6, 59 We will not address such measurements here as their calculation would also require the modeling of the interaction with the conducting substrate (i.e. image charge effects), which goes beyond the scope of the present work.

Table II compares HOMO and LUMO energies for PEN and PFP obtained at DFT and
GW level for an isolated (gas-phase) molecule and for a molecule embedded in the bulk crystal. In order to disentangle the different contributions to the energy levels in the solid state, we report results for different DFT starting points, namely performed either in the gas phase or in the presence of a MM embedding with its proper crystal field. In the following we will adopt the "g" or "e" subscripts to label calculations that are performed for gas-phase  Table I.
We start our analysis by considering results obtained starting from a gas-phase DFT calculation (left columns in Table II). The well-known effect of non-local many-body electronic correlations is the large increase of the HOMO-LUMO gap with respect to the KS value obtained with a functional presenting a small amount of exact exchange (25% in the PBE0 case). 90 As reported in Table II,   within 15%) earlier results from classical polarizable models of atomistic resolution, 49,89,92 and from GW calculations using the polarizable continuum model for embedding. 46 93 We notice that a quantitative account of the gap renormalization in the solid state has been recently achieved also at the DFT level using optimally tuned screened range-separated hybrid functionals. 94 We now turn our attention to the results of calculations that start from a ground state DFT calculations performed for molecules embedded in the MM environment, hence experiencing the microscopic crystal field exerted by the surrounding molecules in the bulk solid.
First, we remark that the HOMO-LUMO gap, either obtained at the KS, GW g or GW e level, is to a very good approximation insensitive to crystal field effects, in contrast to the energies of individual orbitals that are significantly affected (see Table II).
Such a result can be rationalized by considering, to first approximation, the superposition of the quadrupolar fields of MM molecules as a uniform potential acting on the QM region, which implies a rigid shift of all occupied and virtual molecular orbitals with respect to their gas-phase values, with negligible effects on the gap and on neutral optical excitations. This is actually the leading effect in the crystalline materials considered in this work, where the crystal field shifts all energy levels by approximately the same amount ∆ cf ∼ 0.15 and ∼ −0.30 eV in bulk PEN and PFP, respectively. The crystal field shifts ∆ cf in Table II are quantified by the difference between the DFT e |GW e and the DFT g |GW e results. Such level shifts are already present in the starting DFT electronic structure and are then reflected in the subsequent GW calculations, as can be inferred from the data in Table II The inhomogeneity of the crystal field at the atomic scale can also alter the shape and spatial extension of the orbitals, further affecting their energies. Such an effect is illustrated in Figure 2, showing the amplitude of the pentacene HOMO in the gas phase (|φ gas HOM O | 2 ) and its difference with the HOMO in the pentacene crystal. In this case, the effect of neighboring molecules is to stretch the HOMO from the π-conjugated region towards C-H bonds. The relaxation of molecular orbital in the crystal field is found to affect orbital energies by a few tens of meV in crystalline PEN and PFP, although a larger influence is expected in disordered environments or in the case of dipolar molecules. The orbital relaxation in the crystal field does also affect the intermolecular charge transfer couplings and band dispersion, as discussed in Appendix B.
Concerning bulk PEN, our result for the gap in Table II (2.9  Both the polarization and the crystal field effects depend on the shape of the sample, leading to different charged excitations in the bulk and at the crystal surfaces where actually these quantities are experimentally measured. In order to approach photoelectron spectroscopy experiments, we hence explicitly computed quasiparticle excitations for a molecule at the specific crystal surface probed in the experiment we aim at modeling, i.e. the (001) and (100) face for PEN and PFP, respectively. Our results for ionization energies at crystal surfaces are reported in Table III. The magnitude of the gap essentially depends on polarization effects. As we discussed in a very recent paper, the gap is ∼ 0.2 eV larger for a molecule at the surface than in the bulk, as a result of the less effective screening at the interface to vacuum. 42 The 10% decrease of polarization from bulk to surface seems to be characteristic for films of standing elongated molecules as PEN and PFP. Very similar gap values are found at the GW e |DFT g and GW e |DFT e level, confirming that the gap is almost insensitive to the crystal field also at the crystal surface.
HOMO and LUMO levels are instead rigidly shifted by the crystal field. ∆ cf is found to be larger at the surfaces we considered (Table III) than in the bulk (Table II), as already reported on the basis of classical electrostatic modeling. 4,54 We recall that the dependence of ∆ cf on the macroscopic shape of the sample and on the crystal facet in 2D slabs both originate from the conditional convergence of charge-quadrupole interactions, which leads to surface-dependent charged excitations. 3,4 We remark that ∆ cf has opposite signs for PEN and PFP both at the surface and in the bulk. The evolution of the HOMO and LUMO levels in PEN and PFP from the isolated gasphase molecule to crystal surfaces is summarized in Figure 3. We recall that our results for ionization at surfaces do apply also to molecular films on insulating substrate (e.g. SiO x ) with same molecular orientations.
The calculated gas-phase levels (GW g |DFT g with cc-pVQZ basis) are in excellent agree- The effects of the environment are then progressively added, starting from polarization   Figure 5), yet it has been reported that nonlocal correlations lead to GW quasiparticle bands that are up to 20% wider than DFT ones. 12,13,15 The latter effect is expect to lead to a small reduction of the gap by less than 50 meV.
The edges of the densities of states obtained from tight binding calculations, shown as blue areas in Figure 3, define our theoretical estimates of the IP and EA of PEN and PFP crystals, compared to experimental values obtained at the same surfaces in   Exploiting the linear response property of the MM subsystem, the reation field matrix elements corresponding to a given orbitals s i on a given atom can be written as In practice, only the reaction field associated with one s orbital per atom need to be calculated in the limit of large MM clusters, allowing to dramatically speed up the calculation of the extrapolated reaction field matrix. A further reduction of the computational burden could be obtained generalizing the strategy outlined above to s orbitals on different atoms, hence computing fewer s orbitals (e.g. one per molecule) and obtaining the other v reac (s i , β ) elements through the reaction to the dipole field created by appropriate linear combination of charges.
Appendix B: Tight binding band structure from first principles inputs The electronic band structure of PEN and PFP is here described with a tight binding model fully parametrized from first-principles. The model accounts for HOMO and LUMO orbitals of the two molecules in the unit cell and considers a two dimensional lattice corresponding to the herringbone plane, i.e. (001) for PEN and (100) for PFP. Dispersion along the plane normal is neglected, owing to the very small intermolecular hopping terms.
The model is parametrized with orbital site energies from GW e |DFT e calculations at crystal surfaces (Table III) and intermolecular transfer integrals obtained at the DFT level (Table V) Table V and, when employed in the tight binding calculations, they lead to bandwidths in very good agreement with plane-wave results. 59,88 It is interesting to note that the dispersion in PFP is almost one-dimensional, as testified by the flat bands along the X − M and Y − Γ paths in Figure 5(b). This results from the very small transfer integrals for both HOMO and LUMO involving pairs with herringbone "T"-like arrangement (last line in Table V).  (Table III)