Accurate Band Gap Predictions of Semiconductors in the Framework of the Similarity Transformed Equation of Motion Coupled Cluster Theory

In this work, we present a detailed comparison between wave-function-based and particle/hole techniques for the prediction of band gap energies of semiconductors. We focus on the comparison of the back-transformed Pair Natural Orbital Similarity Transformed Equation of Motion Coupled-Cluster (bt-PNO-STEOM-CCSD) method with Time Dependent Density Functional Theory (TD-DFT) and Delta Self Consistent Field/DFT (Δ-SCF/DFT) that are employed to calculate the band gap energies in a test set of organic and inorganic semiconductors. Throughout, we have used cluster models for the calculations that were calibrated by comparing the results of the cluster calculations to periodic DFT calculations with the same functional. These calibrations were run with cluster models of increasing size until the results agreed closely with the periodic calculation. It is demonstrated that bt-PNO-STEOM-CC yields accurate results that are in better than 0.2 eV agreement with the experiment. This holds for both organic and inorganic semiconductors. The efficiency of the employed computational protocols is thoroughly discussed. Overall, we believe that this study is an important contribution that can aid future developments and applications of excited state coupled cluster methods in the field of solid-state chemistry and heterogeneous catalysis.


I. INTRODUCTION
The fundamental energy gap or the band gap (BG) is an important intrinsic physical property of any solid material. In fact, all the materials can be categorized into three general classes depending on the magnitude of the measured band gap energies. Materials with negligible band gap energy (∼0 eV) are characterized as conductors or metallic materials. Materials with very large band gap energies (>15 eV) are said to exhibit insulator behavior, while those having intermediate band gap energies are described as materials showing semiconductor behavior. 1−3 For an N electron neutral system, the fundamental band gap (fundamental BG) is defined as the energy difference between the ionization potential (IP N ) and the electron affinity (EA N ) energies (E) so that BG = IP N − EA N . IP N is the minimum energy that is required to remove one electron from the neutral N electron system (IP N = E N−1 − E N ), while EA N is the minimum energy that is required to remove one electron from the N + 1 electron system (EA N The fundamental band gap is also called transport band gap because it determines the minimum energy necessary to create a charge carrier in the material. 4 In an alternative definition, the band gap refers to the excitation energy of the lowest allowed transition and hence is called optical band gap (optical BG). In organic semiconductors, the energy difference between the fundamental and optical band gaps is about 2−3 eV, and quite commonly photoemission and absorption or reflectance optical spectroscopies can exper-imentally probe both of them. In contrast, for inorganic semiconductors, it is not always possible to clearly distinguish between these two types of energy gaps. 4 The experimentally observed fundamental or optical band gaps are further characterized as direct BGs or indirect BGs depending on whether the experimentally probed excitation process involves vertical (nonadiabatic) or adiabatic excitations. 4 Over the past decades, density functional theory (DFT) in the framework of periodic calculations or nonperiodic embedded cluster calculations has been widely applied in electronic band structure studies in an effort to predict this fundamental property. 1,5−8 In particular, a variety of DFT approaches (one-electron energies or Δ-SCF), linear response DFT approaches like the many-body Green's function approximation (GW), the Bethe−Salpeter equation, as well as time-dependent density functional theory (TD-DFT) provide access to the fundamental and optical BG energy with an accuracy that is strongly functional dependent. 1,5,6,9−11 It has been shown that in periodic DFT calculations the local as well as the semilocal and the generalized gradient approximation (GGA or meta-GGA) density functionals, the most affordable functionals for solids, strongly underestimate the magnitude of the band gap energies in the semiconductor class of chemical systems (more than 1 eV on average). 5,6,8 On the other hand, hybrid functionals (global or screened) are known to perform better. Global hybrid functionals like B3LYP, B3PW91, and PBE0 or screened hybrid functionals (in plane-wave DFT) like HSE are able to lower the mean absolute errors in the band gap calculations to about 0.4 eV. 1,5 It should be noted that these errors are very similar to the errors of the many body perturbation GW approach. 1,5 Meanwhile, a wide variety of specialized functionals have been designed like the SCAN meta-GGA functional that minimizes the self-interaction and improves the efficiency of the band gap calculations for solids toward the accuracy of hybrid functionals and GW methods. 7,8 Similarly, the range separated quantum theory project (CAM-QTP) family of functionals has been proven successful in calculating band gaps and charge transfer excitation energies in benchmark studies. 12,13 Diagrammatically derived functionals based on the random phase approximation (RPA) have been also employed to calculate band gaps, but they are limited to small solid-system sizes due to the high computational cost. 14 Alternatively, electron correlation phenomena can be treated in the framework of ab initio wave-function-based methods. These methods are based on a systematic expansion of the many-particle Schrodinger equation. The most powerful methods are based on coupled-cluster (CC) theory; 15−17 by including higher and higher excitations, CC theory systematically converges to the exact solution of the many-particle Schrodinger equation. Fortunately, the convergence with excitation rank has been proven to be rather fast. Already, CC theory with single and double replacements, relative to a reference determinant (usually, the Hartree−Fock wave function), provides results that are considered the "gold standard" in quantum chemistry. 18,19 Owing to exponential ansatz, the truncated CC methods are strictly size consistent at any level of truncation. 19,20 This guarantees the correct scaling of the CC methods with the system size, thus leading to an accurate description of the relative energies along the potential energy surface or between different electronic states. Hence, CC methods, in principle, are ideal tools to treat chemical phenomena in the field of molecular and solid-state chemistry. 21,22 Extension of the CC methodologies to treat excited state phenomena is conventionally performed in the framework of the linear response (LR) 23,24 or equation of motion (EOM) 18,25−27 theories. A particularly attractive variant of EOM-CCSD has been developed by Nooijen and Bartlett and is called similarity transformed EOM-CCSD (STEOM-CCSD). 26,28 Through a carefully chosen sequence of similarity transformations, one can reach a large manifold of singly excited states by diagonalizing a Hamiltonian that only has the dimension of the hole/particle space (e.g., the same size as the CI-singles matrix). Thus, STEOM-CCSD gives efficient access to many excited states. However, at the same time its intrinsic accuracy surpasses that of EOM-CCSD. 26,28 STEOM-CCSD has shown excellent performance in calculating excited state and magnetic properties on a variety of small and medium size chemical systems. 29 Unfortunately, for many years, the applicability of CC methods was limited to rather small molecules (20−30 atoms max.) due to the steep scaling of the computational time with the system size (i.e., O(N) 6 for CCSD and O(N 7 ) for CCSD(T); N is a measure of system size). Thus, until recently, wave function based correlation methods were not tractable for surface problems or were too demanding to be applied routinely. 30−37 However, several important ground and excited state applications at the level of periodic Monte Carlo coupled cluster (CC) and full configuration interaction (full-CI) exist, 33,34,37 while recently band structure calculations of several semiconductors have been performed at the level of periodic equation of motion based coupled cluster (EOM-CC) theory. 38,39 However, it has been shown that exploiting the short-range character of the dynamical correlation in the framework of local correlation techniques can alter the severe limitations imposed by the steep scaling of the wave-function-based methods with system size and provide access to various surface problems. 40 In this direction, periodic local second-order Møller−Plesset perturbation theory (MP2) has been used to treat a number of surface problems. 21,22,41−44 Similarly, for nonmetallic systems, it has been shown that embedded cluster models that treat the long-range electrostatics and polarization on a molecular mechanics level can lead to an effective reduction of the system size that needs to be treated quantum mechanically. 45−50 The combination of these two approaches has been proven instrumental for performing "gold standard" CCSD(T) level calculations for surface systems. 32,51 It has been shown that, provided that the quantum region is described at the level of the domain-based local pair natural orbital version of the CCSD(T) approach (DLPNO-CCSD-(T)), 52 the adsorption energies for a set of small molecules at the rutile TiO 2 (110) surface can be computed with errors of <0.04 eV (1 kcal/mol) with respect to the available experimental values. 32 Local correlation methods in combination with excited state methods are much less explored. It has been shown that correlated methods for excited states based on either configuration interaction or couple cluster schemes like CIS(D), CC2, or ADC2 can be used to generate excited state PNOs, 53−55 thus leading to state specific PNO−CI and PNO-CC based methods. More recently, we have shown that a ground state PNO scheme can be used in combination with the original equation of motion coupled-cluster techniques EOM-CC and STEOM-CC, 18,19,26 resulting in the backtransformed (bt) bt-PNO-EOM-CCSD and bt-PNO-STEOM-CCSD methods (also referred to as EOM-CC and STEOM-CC), which provide a balanced description for valence, charge transfer (CT), and Rydberg states in optical spectroscopy. 56,57 In particular, bt-PNO-STEOM-CCSD has been able to describe the excited states of a wide variety of chemical systems delivering errors below ∼0.1 eV for valence, Rydberg, and CT states. 56−58 In this work, the performance of the bt-PNO-STEOM-CCSD method to calculate a variety of energy gaps in organic and inorganic semiconductors will be evaluated. For this purpose, a test set of organic and inorganic semiconductors will be employed with band gap energies ranging between 1 and 14 eV. In particular, it will be shown that bt-PNO-STEOM-CCSD alone or when combined with embedded cluster approaches in the case of inorganic semiconductors delivers an efficient computational protocol with strong predictive ability for band gap energies. The performance of bt-PNO-STEOM-CCSD will be evaluated against available experimental data and various Δ-SCF/DFT and TD-DFT methods.
II. THEORETICAL CONSIDERATIONS II.A. Relation between Fundamental and Optical Energy Gaps in Particle/Hole and Wave-Function-Based Methods. As was discussed for an N electron system,

Inorganic Chemistry
where IP N and EA N refer to the vertical ionization potential and electron affinity of the system. In the one-electron picture, these quantities can be replaced by the energies of the highest occupied molecular orbital (HOMO) of the N-electron system and the lowest unoccupied molecular orbital (LUMO) − ε a of the N + 1 electron system 7,8,59,13 such that where uppercase labels refer to many electron quantities throughout, lowercase labels to one-electron quantities. Indices i, j, k, and l refer to doubly occupied orbitals in the reference determinant; a, b, c, and d to virtual orbitals; p, q, r, and s to general orbitals; and σ and θ to spin variables. In practice, orbital energy differences have been proven to be a good approximation in calculating the fundamental energy gaps. At this point, a question arises about how the fundamental energy gap (ΔE F IP−EA ) relates to the optical energy gap (ΔE O S 0 →S 1 ), which is defined as the energy difference between the ground state (S 0 ) energy and the excited state (S 1 ) energy. Let us try to answer this question by recalling in this section the principles of delta self-consistent field/DFT (Δ-SCF/DFT) and time depended density functional theory (TD-DFT).
In Δ-SCF/DFT types of approaches, the ground state energy expressions are evaluated using the Δ-SCF orbitals. 60 Δ-SCF/DFT provides direct access to the fundamental energy gaps. In this approach, the energies of the N, N + 1, and N − 1 electron systems are calculated separately, and they are inserted into the fundamental energy gap relation. Likewise, Δ-SCF/DFT can also be used to calculate optical energy gaps by employing the following relation (also known as purification formula): 60 in which the excited singlet (S 1 ), the triplet (T 0 ), and the broken symmetry (BS) states are given by where ψ and ψ ̅ denote spin-up and spin-down orbitals.
In the conventional TD-DFT approach, one needs to solve the following non-Hermitian eigenvalue problem where Ω is the excitation energy and the vectors X and Y represent the occupied−virtual and the virtual−occupied blocks of the first-order density response matrix. In this scheme, the orbital rotation matrices A and B are written as where (pq|rs) are the two-electron repulsion integrals in the Mulliken notation and (pq|f xc |rs) represents the matrix elements of the exchange−correlation kernel in the adiabatic approximation. The third term in eq 4 and the second term in eq 5 describe the contribution of the HF exchange part of the Kohn−Sham operator. In the limiting cases where c x = 0 or c x = 1, eqs 4 and 5 describe TD-DFT employing pure (e.g., GGA) functionals or TD-HF, respectively.
It can be shown that the optical energy gap ΔE O S 0 →S 1 involving i → a long-range charge transfer (CT), core, or Rydberg excitations in the presence of hybrid functionals (c x > 0) is given by 61,62 as in fact all the exchange type integrals in eqs 4 and 5 practically vanish. In contrast, the (ii|aa) integrals represent Coulomb type integrals that are not small and can have values of several electronvolts. In this case, the orbital energy difference in eq 6 clearly does not approximate the state energy difference to any level of accuracy. By inserting eq 1 into eq 6: it is shown that the excitation energies do not either approximate the IP-EA fundamental energy gap unless a nonhybrid functional is considered. We will come back to that point in the numerical results section, where a numerical comparison will be provided between Delta-SCF approaches, orbital energy differences, and TD-DFT calculations. Alternatively, optical and fundamental energy gaps can be calculated by employing wave-function-based methods in the framework of equation of motion coupled cluster theory. This family of methods provides a natural treatment of ionization potentials, electron affinities, and excitation energies. 18,19,26,56,58,63 It should be noted that, in the Δ-SCF/ DFT approach, different sets of orbitals are used for each state of interest, while as the different states are optimized separately, they can be nonorthonormal. This is however not the case in the case of TD-DFT or the STEOM-CC methods in which only one set of orbitals is used and the computed states are necessarily orthonormal. In the following, we describe briefly the principles of the bt-PNO-STEOM-CC method.
II.B. Principles of bt-PNO-STEOM-CC Method. In the canonical STEOM-CCSD method, 18,19,26 one begins by solving the ground state coupled cluster (CC) Schrodinger equation: H̅ |Φ 0 ⟩ = E 0 |Φ 0 ⟩ by using the similarity transformed Hamiltonian H̅ = e −T̂H e T̂w ith the CC ground state excitation  (8) in which Ŝis the transformation operator in singles and doubles truncation, parametrized in terms of IP and EA solutions.
The most expensive step in the STEOM-CC calculation is the ground state couple cluster iterative process, as well as the iterative EOM-IP and EOM-EA calculations. In the bt-PNO-STEOM-CC scheme, the expensive ground state amplitudes are replaced by those obtained from the linear scaling based pair natural orbitals (DLPNO) CCSD. Hence, in the ground state step, the pair natural orbital machinery (PNO) is used to effectively reduce the virtual excitation subspace in the vicinity of the pair occupied orbitals (9) while in a next step the DLPNO quantities are transformed back into the canonical basis.
These amplitudes are then used to solve the EE, IP, or EA equations in the canonical basis. This leads us to the bt-PNO-EOM-CCSD method. In addition, within this scheme the most expensive EOM-EA and some of the expensive terms of EOM-IP calculations are replaced by the seminumerical chain of spheres (COSX) 64,65 method. Further details regarding the bt-PNO-STEOM and the way it is implemented in the ORCA package have been discussed elsewhere. 56 (11) In this approach, the SOC operator is approximated by the spin−orbit mean field (SOMF) operator, 66 which is an effective one-electron operator that contains one-and twoelectron SOC integrals and also incorporates the spin-other orbit interaction. Similar to the multireference configuration interaction (MRCI) and multireference equation of motion couple cluster (MREOM-CC) approaches, SOC is incorporated into the framework of STEOM-CC and TD-DFT by using the bare untransformed SOC operator in the QDPT step. Hence, in eq 11, H SOC is given by (12) where h SOC (x i ) is the effective mean-field one-electron spin− orbit operator and x i and s(i) refer to the coordinates and spinoperators of electron i, respectively. In the presence of a scalarrelativistic potential, picture change effects on the SOC operator are included. A more consistent approach would be to carry out a many-body similarity transformation of the SOC operator first. However, it has been shown previously that QDPT related methods that make use of the untransformed SOC operator offer a reliable description of SOC effects. 67 −72 By making use of the Wigner−Eckart theorem, eq 11 is reduced to Here, m represents the standard vector operator components.
is a Clebsch−Gordon coefficient that has a single numerical value that is tabulated. It satisfies certain selection Further details on QDPT 66,74,75 and the way it is incorporated in various wave-function-based methods to treat various problems in the field of X-ray spectroscopy and magnetism have been reported elsewhere. 67,71,72,76,77

III. CHOSEN STUDY SET OF SEMICONDUCTORS
The chosen test set of the semiconductors includes the naphtalene (2A), anthracene (3A), tetracene (4A), and pentacene (5A) oligoacene, as well as the α-2T, α-3T, α-4T, α-5T, and α-6T oligothiophene organic semiconductors. In the field of inorganic semiconductors, representative examples from the zinc blende, rock salt, and orthorhombic structural groups are included, namely, ZnS, ZnO, MgS, and MgSe from the zinc blende group; LiF, LiCl, NaCl, LiBr, NaBr, MgO, and MgSe from the rock salt group, as well as GeS from the orthorhombic group. Representative molecular structures from the test set are presented in Figure 1, while a structural analysis is presented in the Supporting Information.
Alkali metal halides are wide band gap direct semiconductors with band gaps >6 eV. In particular, LiF shows the largest known band gap of this class, and it is transparent to short wavelengths of UV radiation, which can be used in UV optics. In contrast, metal-oxide, metal-sulfide, and metalselinide semiconductors show smaller band gaps (<6 eV) than the alkali metal halide semiconductors. Experimental fundamental and optical band gaps are presented in Tables S1 and S2.

IV. EXPERIMENTAL ENERGY GAPS OF THE
SEMICONDUCTOR TEST SET IV.A. Organic Semiconductors. The experimental optical energy gaps of oligoacenes and α-oligothiophenes have been obtained from optical spectroscopy measurements, while the respective fundamental gaps are estimated using experimental values of ionization potentials and electron affinities extracted from photoelectron spectroscopic data. 78−81 IV.B. Inorganic Semiconductors. In the case of inorganic semiconductors, optical and fundamental energy gaps cannot always be extracted separately. In fact, as seen in Table S2 for a given semiconductor, there is a large variety to the reported experimental band gap energies. These energies are closely dependent on the type of the experimental spectroscopic measurement, the experimental conditions, as well as the experimental resolution. Focusing on optical absorption spectroscopic measurements, the experimental band gaps of the test set of inorganic semiconductors are presented in Table  S2. It should be noted that these experimental values have served as reference values in various computational studies. 1,5 Among the various experimentally determined band gap energies, the largest variations are observed for the class of wide band gap semiconductors LiF and MgO. In fact, as seen in Table S3, the variations between the different experimental energy gaps range between 2 and 3 eV. Like in the case of the

V. COMPUTATIONAL DETAILS
All calculations were performed with the Orca 4.1 suite of programs. 87 Optical (ΔE O S 0 →S 1 ) and fundamental (ΔE F IP−EA ) energy gaps were evaluated at the TD-DFT, Δ-SCF/DFT, and bt-PNO-STEOM-CCSD levels of theory. For the ground and excited state energetics, a variety of DFT functionals were employed involving the generalized gradients approximation (GGA) functionals BP86, 88 96,97 and by using deconstructed versions of the correlation-consistent triple and quadrupole-zeta quality (cc-pVXZ, X = T,Q) one-electron basis sets. 98 The calculations were accelerated by employing the resolution of identity approximation (RI) 99 for the Coulomb integrals, while the exchange terms were efficiently computed using the "chain-of-spheres" (COSX) 64,65 approximation by utilizing the def2/J and the cc-pV(T/Q)Z/C auxiliary basis sets, respectively. All the calculations were performed using the secondorder Douglas−Kroll−Hess correction (DKH2) 100,101 to account for the scalar relativistic effects, employing the finite nucleus model. 102 Spin−orbit coupling effects (SOC) were computed in the framework of quasi-degenerate perturbation theory (QDPT). 66,74,75 All the calculated energetics are zero-point energy (ZPE) corrected from respective DFT frequency calculations. For bt-PNO-STEOM-CCSD, an extrapolation to the basis set limit (CBS) was carried out based on the two point coupled-electron pair extrapolation scheme with the all electron cc-pVTZ/cc-pVQZ basis sets. 103 Within this scheme, only the correlation energies of the ground and the excited states are taken into account for the extrapolation of the excited states energies. By contrast, DFT/cc-pVQZ calculations are considered to be at the CBS. Atomic coordinates of all systems were obtained from the American Mineralogist Crystal Structure Database 104 and the Crystallography Open Database. 105,106 VI. EMBEDDED CLUSTER APPROACH The band gap calculations on the inorganic semiconductors were performed in the framework of the embedded cluster calculations. A schematic representation of the employed embedded cluster approach is presented in Figure S1. As shown in Figure 2, all quantum clusters (QCs) were extracted from the crystallographic supercells by either employing a unitcell approach, in which the resulted clusters contain building units that correspond to the crystallographic unit cell, or a ligand-field approach, in which the building units retain the metal−ligand coordination environment. Figure 2 illustrates both approaches for LiCl at the converged cluster size, and it is found that both lead to identical excited state properties. In particular, Figure S4 shows that they converge to the same band gap energy. In order to account for long-range Coulombic forces, point charge fields (PCs) were constructed by including an amount of 5000 to about 10 000 charges. These charges were placed at the appropriate crystal lattice atom positions. In the employed embedding scheme, the positions and magnitudes of the point charges are kept fixed, leading to a nonpolarizable embedding scheme. Between the QC and PC regions, a third boundary region (BR), equipped with capped effective core potentials (cECPs), was introduced, in order to avoid spurious electron leakage from the clusters. In particular, an up to double layer of cECPs, ECP10MDF 107 and ECP2MWB 108 (included in the def2-SD framework), was used to replace the metal and the main group atoms, respectively. The charges that were chosen to equip the cECPs and PCs regions were obtained by (1) imposing cluster neutrality conditions (i.e., constraining the absolute total charge of the system to a value close to zero, as described previously 68 ) and (2) by ensuring uniform charge distribution in both QC, BR, and PC regions. Figure 3 demonstrates both the convergence of the optical band gap of LiCl at the STEOM-CC level with respect to the cluster size and the dependence of the converged result on the input charge of the BC + PC region. The cluster size can then be chosen such that the calculated BG curves converge; in the present case, this happens at n = 4 for all input charges. As for the input charge, the best agreement with experiment as a function of growing cluster size is obtained for q Li = 0.95, which is the very choice that satisfies the criterion of a uniform charge distribution between the QC and BR + PC regions. This value can be obtained by locating the intersection of the reference charge line (qQC = q(BR + PC)) and a line determined using a simple linear-regression between the input charge that equips the BR + PC regions and the average charge distribution of the Li atoms in the QC region (Figure 3b). For this purpose, the electrostatic potential (ESP) charges were used as they are obtained from population analysis after the converged self-consistent field (SCF) iterations. Alternatively, NPA, Hirshfeld, or Loẅdin charge schemes may be used. For ionic systems, Mulliken charges can also be used. However, as this charge scheme is quite vulnerable with respect to the basis set extension, this choice is generally not recommended. The number of point charges, cECPs, as well as the corresponding charges for all clusters is presented in Table S3.   Table 1 and show clearly that both Δ-SCF/DFT and TD-DFT behave similarly in calculating the triplet excitation energies. These energies are in good agreement with the experiment for all the employed functionals. For the singlet excitation energies, TD-DFT predictions are generally higher than Δ-SCF/DFT, while both energies increase in the sequence GGA (BP86, BLYP, PBE), hybrid (B1LYP, B3LYP, BEPW91, PBE0), and double hybrid (B2PLYP) functionals. The two approaches differ by about 0.7 eV for GGA functionals and 0.5 eV for hybrid functionals, while they actually match for double hybrid functionals. One should notice that the orbital energy difference between HOMO and LUMO (Δε HOMO−LUMO ) is an excellent approximation to the average excitation energy between singlet and triplet states av.ΔE Δ-SCF (S 0 →S 1 ,S 0 →T 0 ) only if there is no HF exchange in the functional (BP86, BLYP, PBE). This is due to the fact that in the presence of hybrid or double hybrid functionals the Δε HOMO−LUMO orbital energy difference drastically overestimates the singlet excitation energies and approaches the value of the fundamental energy gap ΔE F IP−EA as the fraction of the HF exchange in the functional increases. In pure HF calculations, Δε HOMO−LUMO becomes higher than 8 eV. In accordance with the discussion in the theory section, the orbital energy differences do not relate directly to either the optical or the fundamental energy gaps.
For the fundamental energy gaps (ΔE F IP−EA ) in general, Δ-SCF/DFT in the presence of hybrid and double hybrid functionals provides good agreement with the experiment. By contrast, the GGA functionals provide rather unbalanced predictions. This is also consistent with previous studies for organic and inorganic semiconductors. 9,78 As expected among all methods tested, the bt-PNO-STEOM-CC provides the best agreement with the available experimental values for singlet and triplet excitation energies as well as the fundamental energy gap. It should be noted that, due to the rigid nature of the oligoacene organic semiconductors, the experimental and calculated energy gaps of the anthracene molecule relate directly to the band gap energies of the actual anthracene organic semiconductor. 111 This holds in general for the majority of the organic semiconductors.

VIII. BAND GAP ENERGIES OF INORGANIC SEMICONDUCTORS: CLUSTER SIZE CONVERGENCE IN THE CASE OF ZNS AND LICL
Entering the field of inorganic semiconductors, we first investigate the convergence of the cluster size with respect to the calculated band gap energies. To this end, we choose to investigate the ZnS and LiCl semiconductors. The experimental optical band gap energies are ΔE O = 3.78 eV and ΔE O = 9.40 eV, respectively (Table S2). For ZnS, PBE periodic band structure calculations have predicted a band gap energy that amounts to 2.60 eV. 1 For LiCl, PBE periodic band structure calculations have predicted a band gap energy that amounts to 7 eV, 1 while B3PW91 periodic band structure calculations have predicted values that range between 8.60 and 8.80 eV. 1,5 In the following, we employ DFT orbital energy differences (HOMO−LUMO gap) as well as TD-DFT and bt-PNO-STEOM-CC state energy difference to calculate the ZnS and LiCl optical energy gaps at various cluster sizes. The results are presented in Figure 4 and Figure S6.
In the case of ZnS, already for cluster [Zn 2 S 7 ] −10 , the PBE calculated orbital energy difference converges to the PBE periodic orbital energy difference. Hence, at the DFT level, the cluster size converges rapidly, and for all intents and purposes, the [Zn 4 S 10 ] −12 cluster can be considered converged. Likewise, in the case of LiCl, the PBE calculated orbital energy difference converges to the PBE periodic orbital energy difference at a cluster size of [Li 8 Cl 8 ] 0 ( Figure 4).
As seen in Figures 4 and S6, TD-DFT shows a similar convergence pattern. Consistent with the results of the previous section, it is observed that, as long as there is no HF exchange in the functional, the orbital energy difference between HOMO and LUMO is almost identical to the energy of the first excited state (ΔE TD-DFT

Inorganic Chemistry
Article that it is a valid choice to identify the converged cluster sizes at the DFT level. As seen in Figure 4, in the case of ZnS and LiCl, the PBE calculated orbital energy differences (obtained by periodic or embedded cluster calculations) or the PBE/TD-DFT state energy differences (obtained by embedded cluster calculations) underestimate the experimental band gap energies by more than 1 and 2 eV, respectively. These errors can be effectively reduced with the use of hybrid functionals. 1,5 In fact, in the case of LiCl, the B3PW91 calculated orbital energy differences (obtained by periodic or embedded cluster calculations) range between 0.8 and 0.6 eV from the experimental values ( Figure S6). Consistent with the analysis described above, the B3PW91 TD-DFT results (obtained by embedded cluster calculations) underestimate the experimental value by about 1.1 eV.
By contrast, in the case of ZnS, the bt-PNO-STEOM-CC results employing the [Zn 4 S 10 ] −12 cluster deviate from the experimental value by less than 0.2 eV. Likewise, in the case of LiCl, the errors from the experimental optical energy gap drop below 0.15 eV only when the bt-PNO-STEOM-CC method is employed together with cluster [Li 8 Cl 8 ] 0 . Compared to the bt-PNO-STEOM-CC results, the PBE or B3PW91 calculated orbital energy differences (obtained by periodic or embedded cluster calculations) range between 0.6 and 0.4 eV, while the TD-DFT/B3PW91 calculated state energies (obtained by embedded cluster calculations) are off by more than 1 eV.
To conclude this part, it has been shown that the combination of the bt-PNO-STEOM-CC method with the embedding cluster approach results in a calculation protocol for energy gaps that shows a relatively fast cluster size convergence while delivering very accurate predictions for the optical band gaps. This is also reflected in Figure S5, in which the orbital energies as well as the excited state energy spectrum are calculated as a function of the system size in the case of LiCl. As seen for small cluster sizes ([LiCl] 0 , [Li 2 Cl 2 ] 0 ), the molecular picture is preserved, while the [Li 8 Cl 8 ] 0 cluster shows a solid behavior which allows for the orbital energies to be assigned either to the occupied or to the unoccupied bands while the excited states start to form an excited state continuum. The final clusters employed in this work are listed in Table S3. In all studied cases, the converged clusters contain less than 10 building units. This relatively fast cluster size convergence may be related to fact that the band gap represents the minimum of the band excitation energies in the solid-state context. Hence, one may reach convergence for this property with cluster size much sooner than convergence of the entire spectrum of excited states. As is seen in Figure S7, in the case of LiCl, the first 10 bt-PNO-STEOM-CC calculated excitation energies are essentially converged at the [Li 8 Cl 8 ] 0 cluster size, while convergence of higher excitations energies (>20) requires much larger cluster sizes.  In the next step of our analysis, the optical band gap energies (Table S1) for the chosen test set of semiconductors are computed in the framework of TD-DFT and bt-PNO-STEOM-CC methods. The results are presented in Table S4 and It should be mentioned that the balance between achievable accuracy and the computational effort is a requirement in any practical study. Thus, focusing on the computational efficiency of the employed methods, the timings (using four CPU cores) in a typical single-point DFT calculation (cc-pVTZ basis) of the Li 8 Cl 8 cluster range between 10 min for B3PW91 to about 1 day for the respective bt-PNO-STEOM-CC calculations. As is seen in Figure S8 within the various calculation steps, the major time is spent in computing the "dressing" integrals, the EOM doubles, and the STEOM-CCSD amplitudes. We have recently extended the EOM-CCSD method in the framework of the domain-based pair natural orbitals (DLPNO). 63 On the basis of the performance of the EOM-DLPNO-CCSD method in calculating electron affinities, it is expected that implementation of the relevant STEOM-DLPNO-CCSD method will be much more efficient than the respective bt-PNO-STEOM-CCSD method. Hence, it is expected that the above-reported timings will be drastically decreased toward values that are more in line with the DFT performance.

X. CONCLUSIONS
In conclusion, a systematic bt-PNO-STEOM-CC and DFT benchmark study on band gap energies in a chosen test set of organic and inorganic semiconductors has been presented. In particular, bt-PNO-STEOM-CC, TD-DFT, and Δ-SCF/DFT were employed to calculate optical and fundamental gaps of oligoacene and oligothiophen organic semiconductors. In addition, it was shown that combining the embedded cluster approach with bt-PNO-STEOM-CC, TD-DFT, or Δ-SCF/ DFT yields a computationally affordable computational protocol that is capable of computing a wide range of band gap energies (1.5−13.5 eV) of inorganic semiconductors. In fact, provided that the charge is uniformly distributed between the QC and the BR+PC regions and that the chosen method operates at the basis set limit the cluster size converges rapidly with respect to the computed band gaps. Among the different methods tested, it was shown that for both organic and inorganic semiconductors, bt-PNO-STEOM-CC provides the best agreement with the available experimental values, resulting in errors that are on average lower than 0.2 eV. Relative to the bt-PNO-STEOM-CC and the experimental values, TD-DFT and Δ-SCF/DFT results are of varying accuracy. They basically follow the expectations derived from Jacob's ladder and are consistent with the periodic DFT results. In this context, double hybrid functionals (e.g., B2PLYP) were found to perform best in a series of functionals consisting of GGA, hybrid and double hybrid functionals. Not surprisingly, it is apparent from our results that increased accuracy comes at the

Inorganic Chemistry
Article price of increased computational cost. For example, simple GGA DFT calculations can be performed on the investigated cluster models in just a few minutes. On the contrary, hybrid and more complicated functionals require hours of computation time. Finally, the bt-PNO-STEOM-CC calculations are more expensive than the respective DFT calculations for these types of systems by about 1 order of magnitude In summary, we have demonstrated, that local correlation methods in conjunction with the embedding approach deliver an excellent quality of band gap energies in a wide variety of semiconductors. In fact, bt-PNO-STEOM-CC is a method with high predictive accuracy; hence, it can be used to predict the band gap energies of unknown materials or to benchmark DFT functionals. Further research is undertaken in our laboratories in an effort to deliver more accurate and more efficient STEOM-CC methods that will be able to treat a large variety of excited state problems in the field of solid-state chemistry and heterogeneous catalysis. Our future development efforts are focused on bringing the STEOM-CC method as close to linear scaling as possible (i.e., in the framework of the domain-based pair natural orbitals, DLPNO) and to improve over the efficiency of the employed embedded cluster approach.

* S Supporting Information
The Supporting Information is available free of charge on the ACS Publications website at DOI: 10

Notes
The authors declare no competing financial interest.

■ ACKNOWLEDGMENTS
A.D., R.I., D.M., and F.N. acknowledge Max Planck Society for financial support. The reviewers of the manuscript are acknowledged for their constructive comments.