High-Mobility Naphthalene Diimide Derivatives Revealed by Raman-Based In Silico Screening

Charge transport in crystalline organic semiconductors (OSCs) is considerably hindered by low-frequency vibrations introducing dynamic disorder in the charge transfer integrals. Recently, we have shown that the contributions of various vibrational modes to the dynamic disorder correlate with their Raman intensities and suggested a Raman-based approach for estimation of the dynamic disorder and search for potentially high-mobility OSCs. In the present paper, we showcase this approach by revealing the highest-mobility OSC(s) in two series of crystalline naphthalene diimide derivatives bearing alkyl or cycloalkyl substituents. In contrast to our previous studies, Raman spectra are not measured, but are instead calculated using periodic DFT. As a result, an OSC with a potentially high charge mobility is revealed in each of the two series, and further mobility calculations corroborate this choice. Namely, for the naphthalene diimide derivatives with butyl and cyclopentyl substituents, the estimated room-temperature isotropic electron mobilities are as high as 6 and 15 cm2 V–1 s–1, respectively, in the latter case even exceeding 20 cm2 V–1 s–1 in a two-dimensional plane. Thus, our results highlight the potential of using the calculated Raman spectra to search for high-mobility crystalline OSCs and reveal two promising OSCs, which were previously overlooked.


Introduction
Crystalline organic semiconductors (OSCs) with efficient charge transport, i.e., with high charge-carrier mobilities µ > 1 cm 2 V -1 s -1 , are a hot topic today, being necessary for the efficient operation of organic electronic devices [1][2][3]. The highest charge mobilities for organic semiconductors were observed in crystalline OSCs composed of π-conjugated molecules with fused heteroaromatic and/or aromatic cores [1,4]; nevertheless, they are low in comparison with the estimations of the physical limit for the intrinsic mobility, i.e., the one in the absence of defects, in OSCs (~70 cm 2 V -1 s -1 ) [5]. Thus, the search for higher-mobility OSCs remains an urgent task for organic electronics. However, reliable measurement of the intrinsic µ in OSCs is not straightforward. For instance, a commonly used method based on fabrication of organic field-effect transistors with single crystals of the studied OSCs as active layers is complicated by the issues with the contacts, the architecture, the dielectric, etc. [6]. When polycrystalline films (which are much easier to process) of the OSCs are used as an active layer, the mobility is considerably (by over an order of magnitude) decreased by grain boundaries and other defects; moreover, relative mobilities for different OSCs do not necessarily represent the trends in the intrinsic mobilities [1,7]. Therefore, theoretical estimation of µ prior to its time-consuming measurements in devices is important for a focused search of high-mobility OSCs. However, though first-principle approaches are generally able to yield relatively accurate mobility predictions, they utilize a number of where e is the elementary charge; a is the lattice constant; k B is the Boltzmann constant; T is the absolute temperature; and the transient localization time τ in roughly corresponds to the characteristic periods of LF vibrations strongly modulating the transfer integrals J. The inverse square of the relative disorder, J/σ J 2 , in Equation (1) represents a typical scaling of the localization length of charge carriers around the room temperature [24][25][26]. In view of the existence of such correlations between µ and σ J , in what follows, we will mainly focus on the search of OSCs with a low dynamic disorder, in pursuit of high mobility. The bottleneck of the σ J and hence the µ calculation is the estimation of the nonlocal electron-phonon interaction (NLEPI), which converts the thermal fluctuations of the structural parameters (intermolecular distances and molecular orientations) into the modulations of the transfer integrals and hence suppresses the mobility [27][28][29][30]. The strength of NLEPI can be quantified by the lattice distortion energy L [31][32][33], which is proportional to σ 2 J in the low-frequency limit, and by the contributions of various vibrational modes to it, L i (for details, see Section 3.3). Very recently, we have suggested a Raman-based approach for estimation of NLEPI [24,31,[34][35][36]. Specifically, we have shown that for LF vibrational modes, L i correlates with the Raman intensity of the mode, I i , yielding an estimate of the total lattice distortion energy: where I(ω) is the spectral Raman intensity (technically, the Raman scattering cross section per one OSC molecule); ω is the vibrational frequency (or the Raman shift); n(ω) is the population of the corresponding vibrational mode; andh is the Planck's constant. The pre-factor C LF = C LF (J, α, E g −hω L , a, . . . ) is expected to depend on the molecular and crystal structure: the molecular polarizability α, the optical or HOMO-LUMO gap E g minus the incident photon energyhω L , the transfer integrals J and lattice constant(s) a, etc. [31]. Equation (2) was used in Ref. [31] mainly to uncover a clear correlation between the experimental Raman intensities of the LF vibrational modes and the calculated L i . In a subsequent work [24], an approach developing Equations (1) and (2) was tested on experimental Raman spectra, which showcased the efficiency of OSC screening based on the experimental Raman + in silico methodology. Note that contemporary Raman spectrometers enable reliable measurements down to~5 cm -1 [37] and even below that, i.e., they can, in principle, detect all the Raman bands of most crystalline OSCs. A theoretical analysis of the physics underlying correlation (2) also revealed a reasonable scaling law of the pre-factor, C LF ∝ (E g −hω L ) 2/ α 2 (which corresponds to the parameter γ = 1/2 in Ref. [24]). Note that in the classical, LF approximation, n(ω) + 1~k B T/hω, and Equation (2) can be simplified: L~C LF I LF /k B T, where I LF = LF I(ω)dω is the total Raman intensity in the LF range. However, both the calculation and measurement of the absolute Raman intensity is not straightforward, with the spectral function I(ω) usually revealed up to a normalization constant (e.g., DFT packages typically use the Raman band with the highest intensity for this constant), thus, it appears advantageous to transform Equation (2) into an expression containing intensity ratios instead. It is convenient, for instance, to use a ratio of the LF Raman intensity to the "high"-frequency one (HF, ω > 200 cm -1 ). The latter is mainly associated with the electron-phonon interaction with intramolecular vibrational degrees of freedom, which is quantified by the reorganization energy, λ. This energy allows for an estimation λ ∼ C HF HF [24,38]; moreover, the physical analysis in Ref. [24] suggests a scaling C HF (α, E g −hω L , a, . . . ) ∝ 1/α 2 (which corresponds to the off-resonant case γ = 1 in Ref. [24]). Now, by dividing Equation (2) over the estimation for λ, one can eliminate the overall normalization of the Raman intensity in favor of a ratio of the LF intensity to the (weighted) HF one: For compounds with a similar chemical structure-including the naphthalene diimide (NDI) derivatives studied in the present paper, see below-reorganization energies are of the same order (see, e.g., Figure 6 in Ref. [34] and Table S7 (Supplementary Information)); the same applies to the energy gaps. Thus, R can be used for estimating L and thus the dynamic disorder amplitude σ J (for details, see Refs. [24,36]): at least, the higher R, the higher L is expected for given λ and E g . The value of R can be estimated both experimentally and computationally. In the present study, we test a Raman-based in silico approach for the search of high-mobility crystalline OSCs, relying upon the Raman intensities calculated by periodic DFT.
A good benchmarking set for such tests is a series of OSCs differing in the alkyl chains end-capping the same π-conjugated unit [39,40]. For instance, NDI derivatives with alkyl and cycloalkyl terminal substituents have been synthesized, with the (polycrystalline) thinfilm OFET charge-carrier mobilities measured for them [2,[41][42][43]. Some of these compounds have shown high (close to record) electron mobilities, for instance, crystalline N,N -dicyclohexyl-naphthalene-1,8;4:5-dicarboximid (NDI-CHex) [41] and N,N -bis(hexyl)naphthalene diimide (NDI-Hex) [42]. However, charge-carrier mobility data for single crystals (i.e., estimates of the intrinsic mobilities) are lacking for most of the NDIs, so intrinsic mobilities for these OSCs are unknown. Thus, application of the Raman-based estimation of the dynamic disorder and calculation of single-crystal charge mobilities for various NDIs can reveal previously overlooked high-mobility OSCs.
Six alkyl-and cycloalkyl-substituted NDIs are considered in this work, including the abovementioned NDI-CHex and NDI-Hex, see Figure 1. We evaluate the Raman spectra of these crystals using periodic DFT, with special attention paid to the accuracy of the calculated wave numbers and Raman intensities for the lowest-frequency Raman-active vibrations.
We find that the lowest R ratios (expected to correspond to the highest charge-carrier mobilities) are observed for NDIs with butyl (NDI-But) and cyclopentyl (NDI-CPen) substituents in the considered family of OSCs. We then proceed with an in-depth investigation of NLEPI in the six crystals and calculate the charge-carrier mobilities (including the anisotropy and temperature dependencies) within the transient localization scenario [11,19]. Note that, in the context of the present work, the role of the mobility calculations is twofold: first, they test the overall performance of the R-based in silico mobility screening approach suggested; second, they yield quantitative mobility predictions for the two shortlisted NDIs-beyond the qualitative conclusion of them being the best ones in the considered series. Namely, our analysis corroborates that these OSCs can show high room-temperature electron mobilities: about 13 cm 2 V -1 s -1 along the high-mobility axis of NDI-But and around 23 cm 2 V -1 s -1 along two axes in NDI-CPen. These values exceed that for the previous record-breaker NDI-CHex, while the one for NDI-CPen even potentially beats the figure for rubrene [44] and, as we further verify, is quite robust against the so-called "killer mode", dropping only to 10-15 cm 2 V -1 s -1 in supercells supporting this mode [45]. Thus, our results both highlight the potential of using the Raman spectrum for screening of high-mobility crystalline OSCs among compounds with similar chemical structure and suggest a potentially record-high-mobility material, NDI-CPen, relating its expected properties to the Raman spectrum. Six alkyl-and cycloalkyl-substituted NDIs are considered in this work, including the abovementioned NDI-CHex and NDI-Hex, see Figure 1. We evaluate the Raman spectra of these crystals using periodic DFT, with special attention paid to the accuracy of the calculated wave numbers and Raman intensities for the lowest-frequency Raman-active vibrations.
We find that the lowest R ratios (expected to correspond to the highest charge-carrier mobilities) are observed for NDIs with butyl (NDI-But) and cyclopentyl (NDI-CPen) substituents in the considered family of OSCs. We then proceed with an in-depth investigation of NLEPI in the six crystals and calculate the charge-carrier mobilities (including the anisotropy and temperature dependencies) within the transient localization scenario [11,19]. Note that, in the context of the present work, the role of the mobility calculations is twofold: first, they test the overall performance of the R-based in silico mobility screening approach suggested; second, they yield quantitative mobility predictions for the two shortlisted NDIs-beyond the qualitative conclusion of them being the best ones in the considered series. Namely, our analysis corroborates that these OSCs can show high room-temperature electron mobilities: about 13 cm 2 V -1 s -1 along the high-mobility axis of NDI-But and around 23 cm 2 V -1 s -1 along two axes in NDI-CPen. These values exceed that for the previous record-breaker NDI-CHex, while the one for NDI-CPen even potentially beats the figure for rubrene [44] and, as we further verify, is quite robust against the socalled "killer mode", dropping only to 10-15 cm 2 V -1 s -1 in supercells supporting this mode [45]. Thus, our results both highlight the potential of using the Raman spectrum for screening of high-mobility crystalline OSCs among compounds with similar chemical structure and suggest a potentially record-high-mobility material, NDI-CPen, relating its expected properties to the Raman spectrum.   [46] are given in parentheses. The red, blue, and gray balls represent the oxygen, nitrogen, and carbon atoms, respectively, while the small white balls indicate the hydrogen atoms.

Raman Spectra of NDI Derivative Crystals, the R Ratios, and Dynamic Disorder
The lowest frequency Raman-active vibrations of single-component organic crystals are often characterized by relatively low Raman-scattering intensities, namely, ~5% of the most intense band [47][48][49][50][51]. Thus, errors in the calculation of wavenumbers (up to 10 cm -1 ) and/or of the Raman intensities (up to two times) of such modes do not play a significant role in the R calculation. In contrast to the general case of organic crystals, however, in OSCs [31,52] and in some multicomponent organic crystals [53,54], the lowest frequency  [46] are given in parentheses. The red, blue, and gray balls represent the oxygen, nitrogen, and carbon atoms, respectively, while the small white balls indicate the hydrogen atoms.

Raman Spectra of NDI Derivative Crystals, the R Ratios, and Dynamic Disorder
The lowest frequency Raman-active vibrations of single-component organic crystals are often characterized by relatively low Raman-scattering intensities, namely,~5% of the most intense band [47][48][49][50][51]. Thus, errors in the calculation of wavenumbers (up to 10 cm -1 ) and/or of the Raman intensities (up to two times) of such modes do not play a significant role in the R calculation. In contrast to the general case of organic crystals, however, in OSCs [31,52] and in some multicomponent organic crystals [53,54], the lowest frequency Raman-active vibration can be very intensive; moreover, it can be located around or below 10 cm -1 [55,56]. For such a mode, errors reaching 10 cm -1 in the calculation of the wavenumber, caused by basis set changes (see Tables S1 and S2 in Ref. [31]) or by the choice of the functional (see Tables S1 and S2), which are observed for the OSCs studied, can radically change the frequency of the mode. Note that 10 cm -1 is also a typical value of the Lorentzian broadening used when comparing vibrational bands with experimental ones, see, e.g., [57]. Moreover, the calculated Raman intensity of the lowestfrequency band can also differ significantly from the experimental value (see Figure S3 in Ref. [58]). Such discrepancies may be due to the limited applicability of the "double harmonic approximation" [59] for the quantitative description of the Raman spectra of organic crystals below 20 cm -1 . For these reasons, in the following, the Raman-active bands below 20 cm -1 are treated with care, and their DFT predictions are not considered as ultimate.
In this regard, especially since LF vibrational modes are concerned, it is worth making some general remarks on the accuracy of DFT calculations of Raman spectra of OSCs, before switching to a full analysis of the six NDI crystals. We start with the calculations of the Raman spectra of crystalline NDI-Hex and NDI-CHex, for which experimental data are available [31], using several functionals and the 6-31G** basis set, in order to choose the one that fits the data best. A comparison between the calculated and experimental spectra is given in Tables S1 and S2. The PBE-D3/6-31G** level gives a better description of the spectra compared to PBEsol/6-31G** and B3LYP/6-31G**, in line with our earlier findings [53]. In particular, we are able to conclude that the wavenumber accuracy for the LF Raman-active bands of the two crystals calculated using PBE-D3 is acceptable: e.g., for the modes above 20 cm -1 , the calculations agree with experiment within 5%.
With the above in mind, we now present the Raman spectra of the six studied OSCs calculated using periodic DFT at the PBE-D3/6-31G** level in Figures 2, 3 and S1 . Their common features are: (i) the most intense Raman band lies in the LF region of the spectrum; (ii) intense Raman bands are practically absent in the 200-500 cm -1 range, in line with Refs. [31,51,54,55] (we use the term "intense band" for vibrations having relative Raman intensities of more than 5%, compared with the most intense band hereinafter normalized to 100%); (iii) the Raman spectra of the crystals under consideration are very similar in the HF range (with the only exception in the range around 3000 cm -1 related to CH vibrations of the substituents), which can be readily explained by the similarity of their conjugated cores, in accord with the literature data (see Figure 6 in Ref. [34]). Thus, the specific features of our crystals manifest themselves in the LF range of the Raman spectrum.
The cycloalkyl NDIs have around 5 intense Raman bands in the LF range (see Figure 2 and Tables S2 and S3). Such spectra, with less or around 5 intense Raman bands in the LF range will be called relatively "sparse" below. The number of intense bands in the LF region of the Raman spectra of the alkyl-substituted NDI crystals is greater than that for the NDIs with cyclic substituents (see Figure 3 and Tables S1-S3), reaching 12 for NDI-Una. This could be related to a softer character of the alkyl chains compared with cycloalkyls; moreover, with the increase in the substituent length, the number of vibrational degrees of freedom also increases.
Finally, the LF/HF ratios, R, obtained from the calculated Raman spectra of NDIs using Equation (3), are summarized in Table 1. In view of the inaccuracies of DFT-calculated Raman wavenumbers and intensities of the lowest-frequency bands discussed above, we additionally estimated the R ratio neglecting the mode below 20 cm -1 (such a mode is present only for NDI-CPen). According to these findings, there are two crystals that are worth particular attention. The first one is NDI-But, which has the smallest R. The second one is NDI-CPen, which has the smallest R if the mode at 8 cm -1 is not taken into account.  Namely, as discussed in Refs. [31,34], the small R values in these crystals should indicate a weak dynamic disorder (and hence, potentially high charge-carrier mobilities) therein. Below, we will analyze how the R values correlate with the dynamic disorder and present the evidence that NDI-But and NDI-CPen are indeed promising for organic electronics. Finally, the LF/HF ratios, R, obtained from the calculated Raman spectra of NDIs using Equation (3), are summarized in Table 1. In view of the inaccuracies of DFT-calculated Raman wavenumbers and intensities of the lowest-frequency bands discussed above, we additionally estimated the R ratio neglecting the mode below 20 cm -1 (such a mode is present only for NDI-CPen). According to these findings, there are two crystals that are worth particular attention. The first one is NDI-But, which has the smallest R. The second one is NDI-CPen, which has the smallest R if the mode at 8 cm -1 is not taken into account.  Namely, as discussed in Refs. [31,34], the small R values in these crystals should indicate a weak dynamic disorder (and hence, potentially high charge-carrier mobilities) therein. Below, we will analyze how the R values correlate with the dynamic disorder and present the evidence that NDI-But and NDI-CPen are indeed promising for organic electronics. In the top panels, contributions from various types of motion (librational and intramolecular vibrational motions of the NDI cores, as well as substituent, "tail" motions) to the vibrational modes are shown (coloring corresponds to Figure 2).
In the NDIs with cycloalkyl substituents, the LF modes with the highest Raman intensities are located below 25 cm -1 and around 100 cm -1 (see Figure 2 and Tables S2 and S3). The atomic displacements for the former modes in these crystals are mainly composed of librations (rotations) of the NDI core around the short axis of the molecule, accompanied by motions of the cycloalkyl substituent (see Table S4; see also Figure 3c in Ref. [31]). The displacements of the atoms of the cyclic substituents are even larger than those of the NDI core atoms in this vibrational mode. In contrast, the higher-frequency modes around 100 cm -1 with dominant Raman intensities correspond mainly to librations of the NDI core around the long axis of the molecule (see Table S4; see also Figure 3e in Ref. [31]). Notably, the atomic displacements of the cyclic substituents are practically negligible in these modes.
In the alkyl-substituted NDIs, the most intense Raman band is located near 100 cm -1 (see Figure 3 and Tables S2 and S3), and the wavenumber of this band in NDI-But, 118 cm -1 , is the highest in the series. Moreover, the wavenumber of the lowest-frequency Ramanactive vibration in NDI-But (63 cm -1 ) is more than twice the corresponding wavenumbers for crystalline NDI-Una and NDI-Hex. Analogously to the cycloalkyl-substituted NDIs, it is possible to single out Raman-intense modes in which the atomic displacements correspond mainly to librations (rotations) of the NDI core around the short axis of the molecule and around its long axis (see Table S4 and also Figure 3d,f in Ref. [31]).
Contributions of librational motions, intramolecular vibrations of the conjugated core, and motions of the (cyclo)alkyl substituents to the vibrational modes are shown in Figures 2 and 3. Unsurprisingly, the lowest-frequency Raman-active modes in the crystals studied are predominantly librational (see also Table S4). Indeed, Γ-point vibrational modes of a crystal with Z red = 1 do not include a translational component, i.e., the centers of mass of the molecules remain at rest. For a perfectly rigid molecule, the residual motion is inevitably rotational, i.e., corresponds to librations around a certain axis. Now, our molecules have a rigid enough condensed aromatic core (which is even more rigid under low-frequency perturbations), attached to flexible (cyclo)alkyl substituents. Thus, it is quite natural that LF Γ-point vibrational modes in Figures 2 and 3 consist of sizeable NDI core librations accompanied by mixed librational/buckling motions of the substituents. Notably, the atomic displacements for librational modes are antisymmetric (odd) with respect to the rotation center, i.e., such modes are Raman-active. The strong Raman signal for these modes (and the strong coupling of these modes to the intermolecular transfer integrals, see below) are expectable, since the corresponding vibrations alter relative orientations of the NDI cores and, hence, the polarizability of the crystal [31]. Strictly speaking, an exception from the above picture is the mode around~125 cm -1 with a considerable intramolecular contribution in all the crystals studied. However, a closer look reveals that this contribution also stems mainly from librations of the central segment of the NDI core around the long axis, which are not synchronized with those of the core edges, leading to its torsional deformation (see Table S4).
Finally, the LF/HF ratios, R, obtained from the calculated Raman spectra of NDIs using Equation (3), are summarized in Table 1. In view of the inaccuracies of DFT-calculated Raman wavenumbers and intensities of the lowest-frequency bands discussed above, we additionally estimated the R ratio neglecting the mode below 20 cm -1 (such a mode is present only for NDI-CPen). According to these findings, there are two crystals that are worth particular attention. The first one is NDI-But, which has the smallest R. The second one is NDI-CPen, which has the smallest R if the mode at 8 cm -1 is not taken into account.  Namely, as discussed in Refs. [31,34], the small R values in these crystals should indicate a weak dynamic disorder (and hence, potentially high charge-carrier mobilities) therein. Below, we will analyze how the R values correlate with the dynamic disorder and present the evidence that NDI-But and NDI-CPen are indeed promising for organic electronics.

Nonlocal Electron-Phonon Interaction
In order to analyze the expected correlation between the computed Raman spectra and the strength of the dynamic disorder (or, equivalently, the lattice distortion energy L, see Equations (4) and (5)), in Figures 2 and 3 we present the contributions of various vibrational modes L i = ∑ mn L i,mn to L and collate them to the corresponding Raman intensities for the crystals studied. As follows from these figures, all the modes with considerable L i values are Raman-active, in line with Ref. [31] and with the above discussion of the Γ-point modes in Section 2.1. Moreover, the contributions of the modes to L correlate with their Raman intensities, also in line with Ref. [31], except for the 97 cm -1 mode in NDI-But having a low Raman intensity but a high L i , and the lowest-frequency modes in NDI-CPen and NDI-CSep, which have very high Raman intensities but weakly contribute to L. As discussed in Section 2.1, the discrepancy between L i and I i for the lowest-frequency vibrational modes could originate from their potentially inaccurate treatment within the used computational methodology.
Being convinced by Figures 2 and 3 that Raman spectrum reveals (to some extent) the contributions of different vibrational modes to dynamic disorder, we turn to checking the correlation between the LF/HF ratio R and the relative dynamic disorder amplitude quantified as σ J /J (for the definition of the latter quantity, see Section 3.3). The data for R, L, and σ J /J, as well as the equilibrium J values, are listed in Table 1; additionally, this table presents the ratios R = λ E g − ω L 2 R/k B T expected to correlate with L (see Equation (3)).
According to the table, the lowest R values are observed for NDI-But and NDI-CPen (for the latter one, neglecting the first Raman-active mode); Table 1 also shows that these crystals demonstrate low disorder amplitudes. Indeed, as follows from the NLEPI details in NDI-But and NDI-CPen given in Tables S5 and S6, the disorder in J is moderate in both crystals: in NDI-But/NDI-CPen, σ J amounts to 0.2J/0.21J in the direction with the largest J (the a/c crystal axis) and to 0.31J/0.15J in the direction with the second-largest J (the b/a axis). An overall correlation between the LF/HF ratios (both R and R) and σ J /J is also vividly shown in Figure S2. One should admit that the correlation between R and the relative disorder (σ J /J) is partially masked by the large transfer integrals, reaching 135 meV in NDI-But and 113 meV in NDI-CPen (cf. the values J = 118 meV in rubrene [60] and J = 69 meV in F 2 -TCNQ [61]-the OSCs showing a record hole mobility and one of the highest electron mobilities reported to date, respectively [1]). However, Table 1 and Figure  S2 show that R also correlates with the absolute disorder amplitude σ J ∝ √ L within each of the two NDI series. Apart from that, another spectral quantity is worth pointing outthe number of intense enough Raman bands in the LF range-and Table 1 demonstrates that this quantity correlates quite well with σ J for the whole set of six NDIs. In any case, based on our analysis, it can be thus expected that (i) NDI-But and NDI-CPen crystals are characterized by potentially high electron mobilities in the considered series of alkyland cycloalkyl-substituted series of OSCs, respectively, and (ii) the quantities R, R (and the number of intense Raman bands) can be used for the screening of high-mobility OSCs. Both these statements are further justified by direct charge transport calculations in the following section.

Charge Transport Characteristics
We finally turn to charge-transport calculations in all the studied NDI derivatives within the transient localization framework (for details, see Section 3.4, and Refs. [11,19,31,62]). As mentioned in Section 1, in the context of the present work, these calculations pursue two goals. First, they are meant to additionally verify our Raman-based screening methodology by comparing the qualitative trends in R and the calculated mobilities (which is potentially important for OSCs in general beyond the two NDI series considered here). Second, after the Raman-based analysis of the previous sections has shortlisted two NDI crystals in these series, NDI-But and NDI-CPen, the calculations are expected to yield quantitative (albeit approximate) predictions of the electron mobilities in them, in cm 2 V -1 s -1 .
Note that for the first of the above goals, it can be instructive to evaluate not necessarily the mobility itself, but also other, more qualitative characteristics relevant to charge transport. For instance, by virtue of the Kubo formula (6), the mobility at a given temperature is proportional to a time integral of the real part of the charge current correlator ĵ a (t)ĵ b (0) , so the profile of the latter can provide signatures of efficient, coherent charge transport in organic crystals [11,19,31]. In particular, the very characteristic correlation decay time τ provides an estimation of the coherence time of charge carriers and, in turn, is able to benchmark their mobility [31]. Thus, along with the calculated mobility values, in Figure 4a-f, we have plotted the current correlators for the six NDI derivatives studied; the calculation has been restricted to the slab plane (xy), as the nonnegligible transfer integrals lie in it, virtually making the transport (at most) two-dimensional (see the discussion of slabs and layers in Section 3.1). One immediately observes that among the three cycloalkyl-substituted NDIs shown in Figure 4a-c, NDI-CPen demonstrates the longest-living correlations; the second-best in this regard is NDI-CHex, and NDI-CSep is the worst, which is well in line with the numbers of intense Raman bands, the relative dynamic disorder magnitudes σ J /J, and the LF/HF ratios R summarized in Table 1. The same applies to NDIs with alkyl substituents (Figure 4d-f): NDI-But is the winner here, characterized with both long-living correlations and their high amplitude along one of the axes (cf. a large transfer integral J 1 and a small J 2 in Table 1); this agrees with the best R ratio, the smallest number of Raman-intense LF bands, and the lowest relative dynamic disorder σ J /J in NDI-But, as compared with NDI-Hex and NDI-Una (see Table 1).
The mobility values resulting from time integration of the correlators also agree with the above semiqualitative analysis based on the correlation amplitudes and decay times. Figure 4a-f gives the eigenvalues µ max , µ min of the 2 × 2 mobility tensor in the slab plane (which are the maximum and the minimum values of the mobility in it, respectively), along with the "isotropic" mobility µ = (µ xx + µ yy )/3 averaged over all directions in 3D space. Here, NDI-CPen, NDI-CHex, and NDI-But demonstrate the highest values, which agree with the first and the last ones being shortlisted from NDIs with cycloalkyl and alkyl substituents, respectively, via a Raman-based analysis. Thus, our results corroborate the suggestion that calculated Raman spectra can be used for assessment of the dynamic disorder and for in silico screening of high-mobility OSC materials.
As for the absolute, quantitative mobility predictions, our estimations are quite close to the µ values available for NDI-CHex and NDI-Hex from both experimental [41] and theoretical [63] studies. For the two champions, NDI-CPen and NDI-But, our calculations demonstrate a pronounced band-like, inverse-power-law temperature dependence of the mobility, µ(T) ∝ T −γ , with NDI-But also having a high mobility anisotropy (see Figure 4g and the inset in it); both these features are attested to be characteristic of high-mobility OSCs with a coherent transport [1]. Notably, for NDI-CPen, which is even poorer in Ramanintense vibrational modes than NDI-(C)Hex, the average electron mobility reaches a recordhigh value of 15 cm 2 V -1 s -1 and the mobility values in the slab plane turn out to exceed 20 cm 2 V -1 s -1 (Figure 4g), a characteristic figure for rubrene single crystals [44]. Finally, since it is known that the effect of the translational (sliding) "killer" mode observed in THz spectra dominates NLEPI and can drastically suppress the mobility in some OSCs [45], we estimated the electron mobility in NDI-CPen, starting from periodic DFT vibrational mode calculations in 1 × 2 × 1 and 2 × 1 × 1 supercells to access such modes. For both supercells, the mobility did drop down, but not drastically, to about 10-15 cm 2 V -1 s -1 , retaining a considerable degree of isotropy, even despite the preferred a or b axes resulting from the choice of the supercell (see Figure S3). Thus, with two directions supporting such mobility values, the predicted average mobility in NDI-CPen is still greater or equal to the figure for rubrene [40], for which the transport is virtually one-dimensional.
In addition, based on our calculations, we are able to conclude that the shorter the (cyclo)alkyl substituent within the series studied, the weaker the dynamic disorder (specifically, the smaller σ J /J) and the larger the charge-carrier mobility, and that cycloalkyl-substituted NDIs show larger µ values than their alkyl-substituted counterparts. Noteworthily, the predicted mobility trends differ from those observed in thin-film OFET measurements in Refs. [42,43]. As we mentioned in Section 1, this can be explained by the fact that in the latter case, charge mobility was limited by grain boundaries and other defects and could be very far from the intrinsic mobility. Thus, we suggest making careful charge mobility measurements in OFETs with NDI-CPen and NDI-But single crystals as the active layers: we anticipate these studies will corroborate high (probably record-breaking) intrinsic electron mobilities in these OSCs. (a-f) and temperature dependences of the mobility in NDI-CPen and NDI-But (g). In panels (a-f), solid/dashed lines show the real/imaginary parts of the correlators; , are the maximum and the minimum mobility values in the slab plane, respectively, and ̅ = ( + )/3 is the mobility averaged over all directions in 3D. These panels also present the characteristic times τ, at which the real parts of the correlators cross zero, as a qualitative estimation of the correlation decay rates. Thick and thin blue lines in the inset in panel (g) represent the mobility anisotropies in the slab planes of NDI-CPen and NDI-But, respectively, with the figures given in cm 2 V -1 s -1 .
As for the absolute, quantitative mobility predictions, our estimations are quite close to the ̅ values available for NDI-CHex and NDI-Hex from both experimental [41] and theoretical [63] studies. For the two champions, NDI-CPen and NDI-But, our calculations demonstrate a pronounced band-like, inverse-power-law temperature dependence of the mobility, ( ) ∝ , with NDI-But also having a high mobility anisotropy (see Figure  4g and the inset in it); both these features are attested to be characteristic of high-mobility OSCs with a coherent transport [1]. Notably, for NDI-CPen, which is even poorer in Raman-intense vibrational modes than NDI-(C)Hex, the average electron mobility reaches a record-high value of 15 cm 2 V -1 s -1 and the mobility values in the slab plane turn out to exceed 20 cm 2 V -1 s -1 (Figure 4g), a characteristic figure for rubrene single crystals [44]. Finally, since it is known that the effect of the translational (sliding) "killer" mode observed in THz spectra dominates NLEPI and can drastically suppress the mobility in some OSCs [45], we estimated the electron mobility in NDI-CPen, starting from periodic DFT vibrational mode calculations in 1 × 2 × 1 and 2 × 1 × 1 supercells to access such modes. For both supercells, the mobility did drop down, but not drastically, to about 10-15 cm 2 V -1 s -1 , retaining a considerable degree of isotropy, even despite the preferred a or b axes resulting from the choice of the supercell (see Figure S3). Thus, with two directions supporting such . Electron transport properties in the slab plane (xy) of NDI derivatives: room-temperature current correlators ĵ x (t)ĵ x (0) , ĵ y (t)ĵ y (0) evaluated within the transient localization scenario (a-f) and temperature dependences of the mobility in NDI-CPen and NDI-But (g). In panels (a-f), solid/dashed lines show the real/imaginary parts of the correlators; µ max,min are the maximum and the minimum mobility values in the slab plane, respectively, and µ = (µ xx + µ yy )/3 is the mobility averaged over all directions in 3D. These panels also present the characteristic times τ, at which the real parts of the correlators cross zero, as a qualitative estimation of the correlation decay rates. Thick and thin blue lines in the inset in panel (g) represent the mobility anisotropies in the slab planes of NDI-CPen and NDI-But, respectively, with the figures given in cm 2 V -1 s -1 .
Summing up, in this paper, we propose a Raman-based in silico screening protocol for the search of high-mobility organic semiconductors, which consists of the following steps. First, an organic compound consisting of a π-conjugated unit with some substituents and characterized by a reasonably good mobility is considered as a starting point. Next, in the structural databases, compounds are searched with a similar chemical structure, that is, a similar π-conjugated unit and various substituents. Preference is given to the crystals that are characterized by one molecule per primitive cell, Z red = 1. Finally, for these crystals, a Raman spectrum is calculated, and for the crystals demonstrating "sparse" LF Raman spectra and/or low LF/HF ratios, the mobility is evaluated using known methods.
Note that the six crystals considered in this work ( Figure 1) were characterized by Z red = 1, unlike the crystals of NDI and naphthodithiophene diimide derivatives with Z red = 2, for which electron mobility has recently been studied, both experimentally and theoretically [66,67].

Periodic DFT Calculations
In the CRYSTAL17 calculations [68], we employed the PBE [69], B3LYP [70,71], and PBEsol [72] functionals with the all-electron Gaussian-type localized orbital basis sets 6-31G**. London dispersion interactions were taken into account by using the semi-empirical D3 scheme [73]. PBE-D3 and B3LYP gave a reasonable description of the LF Raman spectrum of organic crystals in periodic DFT calculations [31,51,53,74,75], while PBEsol was used to study the energetics of molecular crystals [76]. The space groups and the unit cell parameters of the crystals obtained from the X-Ray diffraction experiment were fixed and the structural relaxations were limited to the positional parameters of the atoms [53]. This approximation was widely used in periodic DFT calculations of various properties of organic crystals [77][78][79][80]. Raman intensities were calculated using the "RAMANEXP" keyword of CRYSTAL17 software. The temperature was set to 298 K, the wavelength of the incident laser light 633 nm.
For DFT calculations aimed at estimating the impact of the "killer mode", 2 × 1 × 1 and 1 × 2 × 1 supercells were used, formed by doubling the experimental ones. Further details are given in Section S1.

Evaluation of Electron-Phonon Couplings, Energy Gaps, and Reorganization Energies
The transfer integrals J mn between molecular sites m, n were calculated using the dimer projection method (DIPRO) [81]. Nonlocal electron-phonon couplings ∂J mn /∂Q i , where Q i is the dimensionless vibrational coordinate corresponding to LF mode i, were calculated as the differences between the J mn values in equilibrium and with the atomic displacements corresponding to Q i = 1 along the normal modes from the periodic DFT calculations. The contribution of vibrational mode i with frequency ω i to the lattice distortion energy in the corresponding direction mn, L i,mn [82], was calculated as: and the overall lattice distortion energy in this direction was calculated as L mn = ∑ i L i,mn . The total lattice distortion energy, which includes the contributions from various chargetransfer directions (mn), was obtained as L = ∑ mn L mn = ∑ i,mn L i,mn . The variance of the transfer integrals, which quantifies the dynamic disorder, reads [33]: In line with our previous studies [24,35,36], we estimated the relative standard deviation of the transfer integrals as σ J /J = ∑ mn σ 2 (J mn )/ ∑ mn J 2 mn . For the mobility calculations, we also evaluated the local electron-phonon couplings, ∂ε m /∂Q i , where ε m are the LUMO energies [33]. These couplings were estimated as the differences between the ε m values in equilibrium and with the atomic displacements corresponding to Q i = 1 along the LF vibrational modes. For the evaluation of ε m and J mn , non-periodic (single-molecule and dimer) DFT calculations at the PBE/6-31G** level were performed in the GAMESS package [83,84]. Reorganization energies were calculated using the four-point scheme [33]. Optical gaps were calculated at the B3LYP/6-31G** TDDFT level.

Evaluation of Charge-Carrier Mobility
The charge-carrier mobility tensor was evaluated using the Kubo formula [21]: In this formula, a, b refer to Cartesian axes in the slab plane (see Section 3.1), m, n to molecular sites in a 20 × 20 × 1 "transport supercell" with periodic boundary conditions, and i to vibrational modes on all the sites in this supercell. The operatorĉ n annihilates an electron on site n and x mn is a vector connecting two sites. To yield the charge current correlators ĵ a (t)ĵ b (0) , the time evolution of current operatorsĵ and a further canonicalensemble averaging were computed with respect to the Holstein-Peierls Hamiltonian (see, e.g., Refs. [11,22]): which describes an electron occupying molecular sites in an OSC and interacting with vibrational degrees of freedom. The last, vibrational term, containing phonon creation/annihilation operatorsb † i /b i , was treated quasi-classically, in view of the low frequencies of the vibrational modes studied; such an LF approximation is directly related to the semiclassical transient localization framework [11,19,31,62].
Following a multiscale methodology, the model parameters J mn , ω i , ∂J mn /∂Q i , ∂ε m /∂Q i for a 20 × 20 × 1 supercell were fixed by replication of the transfer integrals, vibrational frequencies, and couplings extracted from periodic DFT calculations of the NDI primitive cells, as well as from single-molecule and dimer DFT calculations (see Sections 3.2 and 3.3). Additionally, to estimate the effect of vibrational modes beyond the Γ point on the mobility in NDI-CPen, periodic DFT calculations in doubled supercells were used instead of those in the primitive cell. Note that LUMO energies ε m are single-molecule quantities and transfer integrals J mn are limited to pairs of neighboring molecules with nonnegligible LUMO overlaps; both are periodic and thus can easily be copied to fill the whole transport supercell. In contrast, calculation of vibrational modes of the whole transport supercell-large enough to contain the partially delocalized charge carrier-using periodic DFT is virtually infeasible, thus, in our transport calculations, LF vibrational modes were treated as dispersionless, localized on specific molecules. For these modes, evaluated using periodic DFT in small cells and further replicated to fill the transport supercell, the couplings ∂J mn /∂Q i , ∂ε m /∂Q i were found.
Temperature dependencies of the mobility were estimated neglecting the changes of the lattice parameters with temperature: the cell was fixed to the experimental one at 293 K. Finally, an exponential cutoff with a characteristic time τ * = /(0.5meV) was added to the current correlators, which is quite standard in transient-localization calculations [11,31]; the particular value of this cutoff, however, has a minor effect on the evaluated mobilities.

Conclusions
Organic semiconductors with "sparse" Raman spectra, i.e., those which consist of a limited number of bands in the LF range, deserve special attention. Namely, as the strength of the dynamic disorder correlates with the Raman LF/HF ratio and with the number of LF Raman bands, these crystals may be characterized by high charge-carrier mobilities. While performing an in silico Raman-based screening of such OSCs, however, one should be careful with the calculated frequencies and the intensities of Raman-active vibrations with wavenumbers below or around 20 cm -1 .
In the present study, using a screening protocol based on the DFT-calculated Raman spectra, we reveal two NDI derivatives with a weak dynamic disorder. Specifically, within the series of alkyl-and cycloalkyl-substituted NDIs, compounds with cyclopentyl and butyl substituents show sparse LF Raman spectra and low LF/HF Raman intensity ratios. Subsequent charge-carrier mobility calculations using transient localization theory indeed predict very high room-temperature electron mobilities in these OSCs: up to 25 and 13 cm 2 V -1 s -1 along certain axes in NDI-CPen and NDI-But, respectively. Thus, our study shows that screening of OSCs based on the calculated Raman spectra can be an efficient approach for searching the materials with a weak dynamic disorder and a high charge mobility, and reveals two promising OSCs for single-crystal n-type organic field effect transistors, one of them characterized by the predicted charge mobility exceeding 20 cm 2 V -1 s -1 . We anticipate that our results will facilitate the search and design of highmobility OSC materials for organic electronics.

Conflicts of Interest:
The authors declare no conflict of interest.