Charge mobility calculation of organic semiconductors without use of experimental single-crystal data

Prediction of material properties of newly designed molecules is a long-term goal in organic electronics. In general, it is a difficult problem, because the material properties are dominated by the unknown packing structure. We present a practical method to obtain charge transport properties of organic single crystals, without use of experimental single-crystal data. As a demonstration, we employ the promising molecule C10–DNBDT. We succeeded in quantitative evaluation of charge mobility of the single crystal using our quantum wave-packet dynamical simulation method. Here, the single-crystal data is computationally obtained by searching possible packing structures from structural formula of the molecule. We increase accuracy in identifying the actual crystal structure from suggested ones by using not only crystal energy but also similarity between calculated and experimental powder X-ray diffraction patterns. The proposed methodology can be a theoretical design technique for efficiently developing new high-performance organic semiconductors, since it can estimate the charge transport properties at early stage in the process of material development.

rank the predicted crystal structures based on their energies, which are calculated using dispersion-inclusive density functional theory (DFT) 8,[24][25][26][27][28][29] . However, the experimentally observed structure does not always show the lowest energy since the structure strongly depends on the crystal-growth conditions, such as temperature, solvents, and film-forming methods. There is also an another practical problem that high computational cost of DFT based methods limits the molecular size and number of crystal structure candidates that can be evaluated.
As for prediction of charge transport properties of organic semiconductors, two conventional models are widely used, namely the incoherent hopping model 30,31 and the coherent band-transport model 32 . However, these models cannot be applied to recent high-mobility organic semiconductors, because experimental studies, such as Hall effect [33][34][35] and electron spin resonance measurements 36,37 , have revealed intermediate transport properties between the hopping and band limits. To overcome this problem, more rigorous theoretical approaches [38][39][40][41][42] including reviews [43][44][45][46] have been recently reported, but the quantitative prediction of the charge transport properties, including temperature dependence, of molecular crystals is still difficult even if accurate static crystal structures are obtained.
Here, we propose a practical method based on low-cost force-field and high-accuracy quantum simulations that can quantitatively estimate the charge transport properties of crystal structure of newly developed molecule from its structural formula using the non-perturbative approach based on quantum wave-packet dynamical simulations coupled with molecular mechanics, which enable us to access the intermediate regime between high-mobility band transport limit and low-mobility hopping transport limit. In the method, the crystal structure data is computationally obtained by searching possible packing structures from the structural formula. Furthermore, to identify the observed crystal structure with practically useful accuracy and low computational cost from a number of computationally suggested polymorphs, we use experimentally measured powder X-ray diffraction (PXRD) data, which has been used for the determination of the crystal structures 47-50 . The PXRD data are easily obtained after early-stage small-amount synthesis. Different from virtual screening approach for a large number of candidates 14,16 , our method combines computational and experimental approaches, and we aim to really create a novel and ideal semiconductor with examinations of a possibility of the synthesized molecule as high-performance semiconductor by the simulations with high-accuracy without a lot of experimental efforts. These theoretical supports can significantly reduce effort and time and provide useful information for the strategic design of single molecules in the development of high-performance semiconductors. Figure 1 shows a schematic diagram of the proposed method for determining crystal structure and carrier transport properties of newly developed molecule. In step (i), we search for energetically favorable conformers based on the structural formula of the target molecule using force field calculations. In step (ii), initial crystal structures are generated using each oriented molecule as the asymmetric unit and common space groups. Then, the initial crystal structures are optimized using molecular mechanics calculations with the force field. The rankings of favorable crystal structures are determined based on assessment value A crystal derived from their calculated crystal energies and similarities between calculated and experimental PXRD patterns. In step (iii), the selected structure according to the ranking is re-optimized using quantum mechanical calculation based on DFT with dispersion correction (DFT-D), and the obtained structure is determined as theoretical crystal structure of Intensity (arb. u) Figure 1. Schematic diagram of the proposed method for determining crystal structure and carrier transport properties of organic semiconductor. (i) Stable conformers of the target are obtained using force field calculations. (ii) Favorable crystal structures are selected from initial candidates based on force field calculations and ranked based on crystal energy and similarity between calculated and experimental PXRD patterns. (iii) Selected crystal structure according to the ranking is re-optimized based on quantum mechanics and theoretical crystal structure of the target is determined. (iv) Hole mobility of the theoretical crystal structure is evaluated using quantum dynamics.
the target. Finally, in step (iv), we evaluate charge transport properties of the theoretical crystal structure using order-N quantum dynamics calculations.
In the present work, we adopt Merck molecular force field 94 (MMFF94) 51 as the force field and DFT-D with the Perdew-Burke-Ernzerhof parametrization of the generalized gradient approximation (GGA-PBE) exchange-correlation functional 52 and the Tkatchenko-Scheffler (TS) scheme 53 (PBE-TS). The similarity between PXRD paterns is calculated by using Gelder's method 54 . Transport calculations are performed using the time-dependent wave packet diffusion (TD-WPD) method 41,55,56 , where the mobility of charge carriers coupled with phonons is evaluated in terms of a time-dependent expression of the Kubo formula using wave packet dynamics. We use 3,11-didecyldinaphtho [2,3- 57 as the practical target for demonstrating our approach and show the procedure for determining the crystal structure and associated carrier transport performance. Since single-crystal structure and charge transport property of C 10 -DNBDT have been well investigated experimentally 58 , C 10 -DNBDT is one of the best choice for validation of our approach.

Results and Discussion
crystal structure determination. In step (i), we found two energetically stable conformers with Ci or C 2 point group symmetries of the C 10 -DNBDT molecule using the MMFF94 potential ( Supplementary Fig. S1). We employed the conformer with the Ci symmetry for the generation of initial crystal structures because known semiconductor molecules usually show a conformation with inversion symmetry in their crystal structure 59,60 . It is considered that for such molecules, this conformation has the advantages of close molecular packing and energetic stability.
In step (ii), 23,808 initial crystal structures were generated using single molecules with different orientations and nine common space groups without inversion symmetry in the Cambridge Structural Database (CSD) 61 . After the optimization of all initial crystal structures using the MMFF94 potential, we found 348 unique structures within 5.0 kcal/mol of the global minimum in crystal energy. Then, in order to identify the observed crystal structure, the unique structures were ranked based on their assessment values A crystal derived from the crystal energy E crystal and S PXRD (Table 1 and Supplementary Table SI). S PXRD shows the similarity between the calculated PXRD pattern of a crystal structure and the experimentally observed PXRD pattern. In the presented method, we assess that the crystal structure with the lowest A crystal value is most similar to the observed crystal structure. According to the ranking based on A crystal , the 1st structure was determined as the observed crystal structure of C 10 -DNBDT in the crystal structure search using MMFF94 potential and was employed for the next step. Here, as references, we show top 4 structures in Table 1 and Supplementary Fig. S3 and discuss them in the supplement. Table 1 shows that we can increase accuracy in identifying the actual crystal structure from computationally suggested structures by using both the crystal energy and PXRD pattern similarity, since the 1st structure has the rank of 25th in the crystal energy and 29th in PXRD pattern similarity. The crystal energy landscapes and assessment value landscapes for C 10 -DNBDT are shown in the supplement to summarize the results of the crystal structure search ( Supplementary Fig. S2).
In step (iii), the 1st structure was re-optimized using DFT-D based on the PBE-TS scheme, and the obtained structure was determined as the theoretical crystal structure of C 10 -DNBDT. In Fig. 2, we superpose the theoretically and experimentally determined crystal and molecular structures. The molecular conformation and packing arrangements in the theoretical crystal structure are in good agreement with those of the experimental crystal structure. These results demonstrate that the proposed method can be used to detect the actual crystal structure of the C 10 -DNBDT molecule.
We next discuss the quantitative difference between the theoretical and experimental crystal structures in terms of Δθ tor , Δθ her , and RMSD 20 in Table 1. θ tor is the torsion angle around bonds between the π-conjugated core and the alkyl group, and θ her is the angle between planes defined by π-conjugated cores of C 10 -DNBDT molecules related by the 2-fold screw axis symmetry in the layer ( Supplementary Fig. S3). RMSD 20 is the root-mean-square deviation in atomic positions, ignoring H atoms, between matching clusters of 20 molecules from the theoretical  and experimental crystal structures 62 . These values were calculated using the software Mercury 63 . The molecular conformation and the herringbone arrangements in the theoretical structure reproduce these features in the experimental structure, with an absolute Δθ tor value of 6.2% and a Δθ her value of −20.6%. The atomic positions in the theoretical structure match those in the experimental structure, with an RMSD 20 of 0.584 Å. In addition, the S PXRD value of the theoretical crystal structure is 0.984 (Table 1). This result indicates that the calculated and experimental PXRD patterns match very well. The effects of the calculation method on the representation accuracy of the experimental structure, which are important for the accuracy of computational determinations for crystal structures and physical properties, are described in the supplement (Supplementary Table SII). Although the MMFF94 and PBE-TS work for estimating C 10 -DNBDT, we agree possibility of existence of more adequate selection of force fields and DFT functionals for organic semiconductors. The selection should be updated as the computational technology advances, from a viewpoint of accuracy and cost. transport property evaluation. In step (iv), we evaluate the hole mobility of a C 10 -DNBDT single crystal using the theoretical crystal structure. A quantum dynamical approach called TD-WPD method 56 is used for this task, where the mobility is evaluated using a time-dependent expression of the Kubo formula. Intermolecular vibrations at around room temperature induce a large dynamic disorder in the transfer integrals and become a major carrier scattering source. Therefore, the effects of the interaction between charge carriers and molecular vibrations on transport properties are introduced as time-dependent transfer integrals. To obtain accurate intermolecular transfer integrals, we evaluate the transfer integrals using DFT-D with maximally localized Wannier functions (MLWFs) 64,65 . The intermolecuar vibrational modes can be obtained with practically useful accuracy based on low-cost force field calculations using the software CONFLEX 66 . The details of the procedure used for the computation of transport properties are described in the methods section and supplement. Before calculating the transport properties, we evaluate the electronic states. A C 10 -DNBDT single crystal is a p-type organic semiconductor, and thus its transport properties are characterized by the band structure originating from the highest occupied molecular orbitals (HOMOs). Figure 3(a-c) show the calculated HOMO band dispersions (left) and intermolecular transfer integrals (right) of the experimental structure obtained at 298 K, those of the theoretical structure, and those of the 1st structure in Table 1, respectively. The intermolecular transfer integrals were obtained as Hamiltonian matrix elements on the MLWF basis set. Because C 10 -DNBDT single crystals with a layered structure are regarded as two-dimensional electron systems, the band dispersions and the transfer integrals on the bc plane are calculated. The transfer integrals of the theoretical structure, shown in Fig. 3(b), are very close to those of the experimental structure, shown in Fig. 3(a), with an accuracy of about 10 meV. The relatively low reproducibility with respect to θ her of the theoretical crystal structure is ascribed to the slight difference in transfer integrals (Table 1). A comparison of Fig. 3(b,c) reveals that the PBE-TS re-optimization of the 1st structure dramatically improved the red-colored transfer integrals, from 9.3 meV to 43.8 meV, making them similar to those of the experimental structure (56.9 meV). This improvement is due to the lattice constant c of 7.048 Å for the 1st structure being overestimated compared to that of 6.112 Å for the experimental structure.
Using the obtained transfer integrals, we estimate the intrinsic mobility along the column direction (c-axis) of a C 10 -DNBDT single crystal. First, we compare the estimated mobility using the theoretical crystal structure with the experimentally observed mobility of a single-crystal field-effect transistor (FET) 58 , shown by red circles and triangles in Fig. 4, respectively. The estimated mobilities, which are 15.7 cm 2 V −1 s −1 at 300 K and 23.6 cm 2 V −1 s −1 at 220 K, show good agreement with the FET mobilities, which are 17 cm 2 V −1 s −1 at 290 K and 20 cm 2 V −1 s −1 at 220 K. The carrier transport exhibits band-like transport, with an increase in mobility with decreasing temperature. This implies that the wavefunctions of holes in a C 10 -DNBDT single crystal spread over the molecules, leading to coherent band-like transport. The discrepancy between the estimated mobilities and the FET mobilities weakly increases with decreasing temperature. This is because actual organic single crystals inevitably have static disorder, caused by defects and impurities, which becomes another carrier scattering source different from molecular vibrations. In contrast, to estimate the intrinsic transport properties, our simulations take only the interaction between electrons and molecular vibrations into account.  Table 1  www.nature.com/scientificreports www.nature.com/scientificreports/ To remove the static disorder effects from the experimentally observed mobility, we calculated the mobility using the experimental structure. The temperature dependence is shown as black circles in Fig. 4. The calculation results imply that the experimentally observed FET mobility (shown as triangles) has the potential to improve up to the values shown as black circles when static disorder is reduced. Furthermore, the slight difference between    www.nature.com/scientificreports www.nature.com/scientificreports/ the red and black circles confirms the sufficient accuracy of the theoretical crystal structure in terms of transport properties. The mobility calculated from the 1st structure (blue circles) becomes much lower than the estimated mobility of the theoretical crystal structure due to the low reproducibility of herringbone arrangements in the force field calculation (Table 1). Therefore, we can conclude that PBE-TS re-optimization is required to increase the accuracy of the theoretically estimated mobility.
computational and experimental cost. The proposed methodology can detect crystal structure and intrinsic transport properties from structural formula and measured PXRD data of a target molecule usually within 1-3 months. The calculation times for steps (i) to (iv) in Fig. 1 for the target molecule considered here are as follows. It took 22 days to find energetically favorable crystal structures using parallel computing with 200 cores. The crystal structure search scales almost linearly with the number of CPU cores 67 ; for example, when using 400 cores, the calculation time can be shortened by about half to about 10 days. The experiments to measure the PXRD pattern can be progressed simultaneously with the theoretical crystal structure search. In many cases, the synthesis, creation of powder crystal and measurement of PXRD data for the target molecule require 1-30, 1-90, and ≤1 days, respectively. Re-optimization using DFT-D required 15 days using parallel computing with 8 cores. The calculation time required to evaluate the transport properties at a given temperature was about 4.5 days, including the generation of material parameters, such as intermolecular transfer integrals and coupling constants between electron and molecular vibrations, using parallel computing with 20-80 cores. Detailed information about the computation times and resources for each calculation procedure is given in the supplement (Supplementary Table SIII). In general, much more time and effort are required to develop practically useful semiconductor material using only experiments because intrinsic charge transport properties of newly developed material are not clear until a experimental measurement for the single crystal is performed and, furthermore, the observed performance is not always high. By applying the proposed methodology to the development of new semiconductor molecules, candidate molecules that theoretically have high performance can be efficiently found at early stage in the process of material development, significantly reducing the experimental effort and accelerating the development process. In order to establish more high-efficient methodology, we should further improve the computational cost by applying high-throughput simulation techniques to our methodology.

Summary
We proposed a method for determining the crystal structure of organic molecule and associated transport properties without use of experimental single-crystal structure data. First, we find energetically stable conformers from the structural formula of the target molecule and obtain favorable crystal structures from the initial candidates based on force field calculations. Second, based on the experimentally observed PXRD pattern and calculated crystal energy, possible crystal structure is selected. Third, theoretical crystal structure of the target is determined by re-optimization using DFT-D based on quantum mechanics. Finally, we evaluate the hole mobility for the theoretical crystal structure via quantum dynamics using the TD-WPD method, where the transport properties can be evaluated by the Kubo formula using the wave packet dynamics of charge carriers interacting with molecular vibrations. As a demonstration, the method was applied to the C 10 -DNBDT molecule, and the crystal structure and its transport properties were determined. Crystal structure candidates were narrowed down from 23,808 initial candidates to 348 using molecular mechanics, and to 1 using PXRD data and quantum mechanics calculation. The theoretical crystal structure and transport properties show good agreement with experimental results. Although we should further benchmark our methodology by using other organic semiconductors to elucidate a usefulness and transferable ability, we showed that our methodology can detect the crystal structure and intrinsic transport properties of practical organic semiconductor molecule from structural formula and measured PXRD data on a realistic timescale with useful accuracy. The strategy in our methodology is available for other targets with structural formula and PXRD data and is expected to be useful for developing new high-performance semiconductor materials.

Methods
procedures of crystal structure determination. The crystal structure determination of an organic molecule is performed in three steps, corresponding to steps (i)-(iii) in Fig. 1. In the present work, the target molecule for practical demonstration was C 10 -DNBDT. The calculations for finding and assessing crystal structures in steps (i) and (ii) are performed using our original method 8,67 and MMFF94 51 implemented in the software CONFLEX 66 . We employ DFT-D with the GGA-PBE exchange-correlation functional 52 , on-the-fly ultrasoft pseudopotential, and TS scheme 53 . Refinement of the crystal structure using DFT-D based on the PBE-TS scheme in step (iii) is executed using the software CASTEP 68 in Materials Studio 2017 69 . The MMFF94 and PBE-TS have been used for our previous studies and have provided a good accuracy for structure reproducibility and relative stability of crystals 8,67,70 .
In step (i), to detect a stable conformation of alkyl groups with an all-trans conformation relative to the π-conjugated core in C 10 -DNBDT, we perform a conformational search by changing the torsion angles around bonds between the π-conjugated core and the alkyl groups in steps of 10 degrees. Here, we employ the all-trans as initial conformation of alkyl groups because the conformation is thought to provide a high density packing as observed in known crystal structures 59,60 .
In step (ii), one conformer obtained from the previous step is selected and rotated around the x, y, and z axes in steps of 30 degrees. The number of molecules in the asymmetric unit is set to one, and nine common space groups without inversion symmetry are employed, namely P1, P2 1 , C2, Pc, Cc, P2 1 2 1 2, P2 1 2 1 2 1 , Pca2 1 , and Pna2 1 , considering the inversion symmetry of C 10 -DNBDT (Supplementary Fig. S1) and space group frequency ranking using the CSD 61 . The initial crystal structures are generated using each oriented molecule as the asymmetric unit and the common space groups as well as some of molecular positions and different lattice constant parameters. (2020) 10:2524 | https://doi.org/10.1038/s41598-020-59238-2 www.nature.com/scientificreports www.nature.com/scientificreports/ Then, we optimize intramolecular geometry, molecular orientation and translation, and unit cell dimensions in each initial crystal structure by minimizing crystal energy E crystal under the specified space group symmetry. Thus, in the optimization, we employ a fully flexible molecular model and the initial molecular conformation including the alkyl groups changes to adequate conformation for packing arrangements due to intermolecular interactions. For all optimized crystal structures, it is confirmed that they have no imaginary frequencies by performing a normal mode analysis in the atomic coordinates under the specified space group symmetry.
Next, unique crystal structures within 5.0 kcal/mol of the global minimum in crystal energy are extracted from all optimized structures because the observed crystal structures were nearly always found within 5.0 kcal/ mol in previous studies 8,70,71 . Here, we introduce an assessment value calculated as A crystal = ΔE crystal + αΔS PXRD , where ΔE crystal is the crystal energy difference from the global minimum, and ΔS PXRD is obtained by subtracting S PXRD from 1.0. S PXRD represents the similarity between the calculated and experimental PXRD patterns and is calculated by using Gelder's method 54 . With α set to 25 kcal/mol, the discrepancy between the calculated and observed structures defined by a ΔS PXRD value of 0.2 is estimated as an energetic penalty of 5.0 kcal/mol in the assessment of possible crystal structures. The parameter α should be adjusted if we change the force field because the effective range of relative energy between polymorphs depends on the force field. The rankings of the unique crystal structures are decided based on their assessment values, A crystal ; the top four structures for the target molecule in this study are shown in Table 1. The details of the procedure in step (ii) are described in the supplement.
In step (iii), the crystal structure, which is selected according to the ranking, is re-optimized using DFT-D with the PBE-TS scheme. For the re-optimization of 1st structure in this work, the cutoff energy for the plane wave and the total energy convergence tolerance are set to 450 eV and 0.5 × 10 −5 eV/atom, respectively, and the Brillouin zone integration is performed with a 1 × 3 × 4 k-point set.
calculations of carrier mobility. We evaluate the intrinsic hole mobility of a C 10 -DNBDT single crystal using our order-N quantum dynamics simulation technique called TD-WPD 56 . The charge transport properties are characterized by the intermolecular transfer integrals, the molecular vibrations, and the coupling between holes and molecular vibrations. The transfer integrals between molecules are computed using the high-accuracy Wannier method 64,65 based on DFT-D 72 . We employ the 1 × 2 × 3 super-cell structure. The cutoff energies for the plane wave and charge density are 80 and 800 Ry, respectively. Brillouin zone integration is performed at the Γ point. The normal modes of molecular vibrations are obtained using force field calculations with low computational cost using the software CONFLEX 66 . The magnitudes of coupling constants are evaluated from the change in transfer integrals due to atomic displacement along the normal mode. We employ a monolayer consisting of 200 × 200 unit cells to obtain the trajectory of wave packets. The wave packet dynamics are computed up to 2 ps with a time step of 0.5 fs. The details of the procedure for the computation are described in the supplement.
Acquisition of experimental reference data for the crystal structure determination. First, we describe the acquisition of experimental PXRD data, which are used for the crystal structure determination. In the preparation of powder samples, oDCB solution of C 10 -DNBDT (2 mg mL −1 ) was slowly cooled from 110 °C to room temperature for 10 h. The resulting crystals were collected by filtration and dried in vacuo at 140 °C. PXRD measurement of C 10 -DNBDT was performed using synchrotron light source with a wavelength of 1.08 Å at BL44B2 of SPring-8 (Proposal 20160029). Microcrystalline powder samples in capillary tube with a diameter of 0.3 mm were used for transmission method with a exposure time of 210 seconds. Diffraction data were collected on an imaging plate.
Second, we describe the acquisition of experimental crystal structure data of C 10 -DNBDT. In order to validate the theoretical crystal structure determined by the proposed methodology and the measured PXRD data, we perform single-crystal X-ray diffraction (XRD) measurement and determine the crystal structure. Single crystals were obtained by means of recrystallization from certain organic solvents. Those C 10 -DNBDT (10 mg mL −1 ) were grown via the gradual cooling of o-dichlorobenzene (oDCB) solution from 120 °C to room temperature for 24 h. The prepared crystals were collected and dried on filter paper. The XRD data for the single crystals of C 10 -DNBDT were collected on a Rigaku XtaLAB Synergy imaging plate diffractometer with Cu K α radiation.
Device fabrication and characterization of c 10 -DnBDt fets. A device was fabricated on a polyethylene naphthalate (PEN) substrate. The PEN (Teijin, Q65HA-125) film was prebaked at 150 °C for three hours in a vacuum oven to reduce thermal expansion, and then cleaned with acetone and 2-propanol in an ultrasonic bath. A 50 nm thick layer of gold as the gate electrode was deposited using thermal evaporation through a shadow mask. As a gate dielectric, a fluorinated polymeric insulator (EPRIMA AL-X601, Asahi Glass Co.) was used. A single crystal of C 10 -DNBDT was then grown on top of the gate dielectric layer from purified 3-chlorothiophene (Pi-Crystal Inc.) with a concentration of 0.025 wt% via continuous edge casting at 110 °C. A gold source and drain electrodes were thermally evaporated through a shadow mask. Finally, laser etching was performed (V-Technology Co., Ltd., Callisto (266 nm)) for the C 10 -DNBDT layer to form a precise Hall bar architecture. The channel length (L) and width (W) were 200 μm and 50 μm, respectively. Four probes were mounted between the source and drain electrodes, where the distance between two longitudinal probes along the channel length (defined as L*) was designed to be 75 μm. The temperature T dependence of the four-terminal mobility measurements was evaluated in a He gas-exchange cryostat. The transfer characteristics with V D = −2 V were continuously recorded by sweeping V G during the slow ramping of T down and up, where the T ramping rate (0.2 K min −1 ) was much slower than the rate for V G scans (approximately 1 V s −1 ).