Accelerated Discovery of High-Refractive-Index Polyimides via First-Principles Molecular Modeling, Virtual High-Throughput Screening, and Data Mining
We present a high-throughput computational study to identify novel polyimides (PIs) with exceptional refractive index (RI) values for use as optic or optoelectronic materials. Our study utilizes an RI prediction protocol based on a combination of first-principles and data modeling developed in previous work, which we employ on a large-scale PI candidate library generated with the ChemLG code. We deploy the virtual screening software ChemHTPS to automate the assessment of this extensive pool of PI structures in order to determine the performance potential of each candidate. This rapid and efficient approach yields a number of highly promising leads compounds. Using the data mining and machine learning program package ChemML, we analyze the top candidates with respect to prevalent structural features and feature combinations that distinguish them from less promising ones. In particular, we explore the utility of various strategies that introduce highly polarizable moieties into the PI backbone to increase its RI yield. The derived insights provide a foundation for rational and targeted design that goes beyond traditional trial-and-error searches.
I. INTRODUCTION<!>II. BACKGROUND, METHODS, AND COMPUTATIONAL DETAILS<!>B. Refractive Index Modeling Protocol<!>C. Molecular Design Space<!>D. Data Mining and Pattern Analysis<!>III. RESULTS AND DISCUSSION<!>IV. CONCLUSIONS<!>SUPPLEMENTARY MATERIAL<!>COMPETING FINANCIAL INTERESTS
<p>Polyimides (PIs), shown in Fig. 1, are synthesized by polycondensation of R 1 -containing dianhydride with R 2based diamine or diisocyanate [1]. They are an appealing class of organic materials due to their exceptional thermal stability and easy processability [2][3][4]. These properties are complemented by mechanical stability, flexibility, light weight, low cost, as well as flame and radiation resistance, and they thus hold much promise for a range of applications [5,6]. However, their generally low index of refraction (RI) undermines their utility for use in many optic and optoelectronic devices [7][8][9][10][11], such as (image) sensors [12,13], displays [14], and light sources (including organic lightemitting diodes) [15], in which organic materials can be deployed in situ as microlenses [16], waveguides [17], microresonators [18], interferometers [19], anti-reflective coatings [20], optical adhesives [21], and substrates [22]. Most of these applications demand large RI values, often upwards of 1.7 or 1.8. Typical carbon-based polymers -PIs included -only exhibit values in the range of 1.3-1.5 [11]. This situation provides a strong incentive to pursue novel high-RI PIs that are suitable for the aforementioned applications.</p><p>As the properties of organic materials can be tailored and tuned by controlling their molecular structure, they are a prime target for combinatorial search or rational design attempts [23][24][25]. The addition of highly polarizable moieties that do not have extensive π-electron conjugation can increase the RI of polymers without negatively affecting other optical properties. (Significant π-conjugation can lead to large optical dispersion and birefringence [26], as well as poor transparency and coloration.) For example, the incorporation of small aromatic rings, halogens, metals, and in particular sulfur have shown promise for this purpose [3,[27][28][29][30][31]. In 2007, Ueda et al. developed PI films with relatively high RI values of up to 1.75, but they, unfortunately, exhibited large birefringence [32]. In recent years, these findings have been improved upon by increasing the sulfur content of the PIs, yielding RI values of up to 1.76 and smaller birefringence [33].</p><p>In this paper, we present a computational and datadriven approach to study the RI values of PIs and rapidly identify promising lead compounds. We investigate different ways that introduce highly polarizable moieties into the PI framework in order to overcome the technical limits of existing compounds. Given the encouraging results from the earlier work discussed above, we put a focus on incorporating sulfur into the PIs, and in doing so, create a new class of high-RI polymers. In Sec. II, we detail the methods employed in this work, i.e., we provide background on our data-driven in silico approach (Sec. II A), introduce our RI modeling protocol (Sec. II B), discuss the molecular design space we consider (Sec. II C), and outline our data mining and analysis techniques (Sec. II D). Sec. III presents and discusses the results of our study, and our findings are summarized in Sec. IV.</p><!><p>A. Data-Driven In Silico Approach</p><p>The development of new materials such as PIs has traditionally been an experimentally-driven, trial-and-error process, guided by experience, intuition, and conceptual insights. This approach is, however, often costly, slow, biased towards certain domains of chemical space, and limited to relatively small-scale studies, which may easily miss promising leads (both on individual compounds as well as compound classes).</p><p>The study at hand instead embraces a data-driven in silico research paradigm that has gained considerable interest in the past few years [34,35] for its promise to address the inherent complexities of structure-property relationships and the vastness of chemical space more efficiently (see, e.g., Refs. [36][37][38][39][40][41][42][43]). Our work brings together molecular modeling, high-throughput computational screening, and machine learning as well as a corresponding software ecosystem to support data-driven discovery and rational design [44,45]. It has its greatest utility as part of integrated research pipelines with experimentalist partners, where it provides guidance for experimental efforts and mitigates some of their beforementioned shortcomings.</p><p>Our research approach and its rationale can be summarized as following: We first establish a computational modeling protocol that can make sufficiently accurate and fast predictions for the target property in the compound class(es) of interest -in the study at hand for the RI values of PIs (cf. Sec. II B). We then create a largescale virtual library of candidates within that compound class (cf. Sec. II C), on which we apply this modeling protocol to evaluate the performance potential of each candidate compound. In addition to obtaining information about each individual compound, we can mine the screening results in their entirety to reveal underlying structure property-relationships (cf. Sec. II D). Our approach can thus identify both specific lead compounds as well as high-value molecular patterns and features for the de novo design of tailored PIs or the creation of additional screening libraries.</p><!><p>According to the Lorentz-Lorenz equation, the RI (n r ) is a function of the polarizability (α) and the number density (N ), i.e., n r = 1 + 2αN/3ǫ 0 1 − αN/3ǫ 0 .</p><p>In previous work [46,47], we developed a modeling protocol that allows us to accurately and efficiently predict the RI values of polymers within the Lorentz-Lorenz equation. This protocol is based on a combination of first-principles quantum chemistry calculations and data modeling. We compute (static) polarizability values α using the coupled-perturbed self-consistent field equations within the Kohn-Sham density functional theory (DFT) framework with the PBE0 hybrid functional [48] and the double-ζ-quality def2-SVP basis set by the Karlsruhe group [49]. We employ an all-electron, closed-shell approach and include Grimme's D3 correction [50] to account for dispersion interactions. The polarizability calculations are performed on geometries optimized using the universal force field (UFF) [51] following a 3D conformational screening as implemented in the OpenBabel software [52]. All DFT calculations are carried out using the ORCA 3.0.2 quantum chemistry package [53]. We compute the number density N via the van der Waals volume V vdW and packing factor K p of the amorphous bulk polymer as</p><p>We obtain V vdW using Slonimskii's method detailed in Ref. [54] and K p using the support vector regression (SVR) machine learning model introduced in Ref. [46] (except during the prescreening of individual residues R 1 and R 2 discussed in the following section, for which we simply employ a constant K p of 0.75 that is typically found in PIs [55]). We obtain the asymptotic limit for the polymers via a linear extrapolation scheme as outlined in Ref. [46]. (In expoloratory calculations, we already observe near-perfect extensivity based on the monomer units, which we explain with the large size of the PI monomers (∼ 150 atoms) and the finite correlation length exhibited in them. There is thus no need to perform expensive oligomer-sequence calculations.) We previously demonstrated on a set of 112 non-conjugated polymers with experimentally known RI values that this protocol can make rapid and accurate predictions [46,47], and thus enable the high-throughput screening study at hand. We execute the latter using our automated virtual screening code ChemHTPS [56], which creates inputs, executes and monitors the calculations, parses and assesses the results, and which extracts and post-processes the information of interest. Where applicable, we report the mean absolute error (MAE) and the root mean squared error (RMSE) as well as their percentage errors (MAPE and RMSPE, respectively) compared to available experimental data.</p><!><p>In Fig. 2, we show the RI heat map in the α/N parameter space and mark the positions of the 112 polymers we previously studied [46,47]. The figure emphasizes the relative importance of the counteracting parameters and their feasibility in real-world compounds. In order to design high-RI systems, we in principle have to pursue compounds that simultaneously feature both large polarizability and number density values. The positioning of the known compounds with high-RI leans more towards high polarizability and relatively small number density rather than vice versa. One path to narrowing down compound space is to maximize the polarizability for systems from the same compound family, for which the number density is similar. Given that this is known to be the case for PIs, we pursue this path and focus our search on highly polarizable PI candidates [55]. For this, we create a PI library based on 29 building blocks and bonding rules shown in Fig. 3, which we select following empirical expertise with regard to promise and synthetic feasibility. The building blocks can be divided into linkers (B1-B6) and (hetero-)aromatic moieties (B7-B29). Using our ChemLG molecular library generator code [57] with a systematic combinatorial linking scheme, we initially generate 38,619 R 1 structures and 171,172 for R 2 . Combining all R 1 and R 2 to form PIs would lead to a total of 6.6 billion compounds. To restrict the search space to a more manageable number of candidates, we select only the most promising R 1 and R 2 groups (which ultimately distinguish the PI candidates) based on the computed RI values of the individual R 1 and R 2 . After prescreening the residues, we develop PI candidates by introducing the top 100 R 1 and 100 R 2 structures into the PI core motif as shown in Fig. 1. We subsequently screen the 10,000 resulting PI candidates. ChemLG keeps a record of the list of building blocks and connections that are used to create a particular compound, and we use this information for the analysis of structure-property relationships.</p><p>FIG. 3. Molecular building blocks used to create the residues R1 and R2 of the PI screening library studied in this work. The "R" on the building blocks denote allowed sites for linking. The building blocks marked in green (B1-B6) are linkers and the blue ones (B7-B29) are (hetero-)aromatic moieties.</p><!><p>In addition to identifying the best candidates from our high-throughput screening, i.e., those with the highest predicted RI values, we further analyze the collected data to develop a better understanding of the structural patterns that lead to high-RI PIs. For instance, we conduct a hypergeometric distribution analysis, in which we compute the Z-scores (Z i ) of each building block i used in the creation of the screening library as</p><p>, where M is the total number of compounds in the entire library, m is the subset of compounds that is considered (e.g., the top compounds), K i is the number of occurrences of building block i in M molecules and k i its occurrences in the subset of m molecules. A large Z-score thus indicates a statistical overexpression of the associated building block in the high-RI candidates relative to the overall screening library. By applying the Zscore analysis, we can thus identify the most important building blocks and the degree to which they correlate with large RI values. We perform a similar analysis of building block combinations to reveal synergistic effects similar to those known from push-pull or donor-acceptor copolymers. In addition, we present an analysis of Zscore trends for each building block in ranked candidate subsets as well as of the average RI values of the candidates derived from each building block. All data mining work was conducted using our ChemML code [58].</p><!><p>We previously tested the RI modeling protocol described in Sec. II B for the before-mentioned set of 112 non-conjugated polymers and established its predictive performance. However, this benchmark set did not contain any PIs. To prove our protocol's validity for this particular compound class, we perform additional calculations on ten PIs with experimentally known RI values. We compare the experimental and computational results in Tab. I. The MAE (MAPE) and RMSE (RMSPE) values for this comparison are 0.021 (1.2%) and 0.025 (1.4%), respectively, which suggests that our protocol can accurately predict the RI values of PIs as well.</p><p>The results of the R 1 and R 2 residue prescreening are summarized in Fig. 4 (a) and (b), respectively, and the results of the subsequent PI screening are given in Fig. 4 (c). The plots show histograms of the computed RI results. Most of the candidates for R 1 and R 2 have RI values between 1.5 and 1.7 with an average of (a) 1.600, and (b) 1.627, respectively. However, the PI candidates derived from the top R 1 and R 2 residues have significantly larger RI values with an average of 1.843. We note that either of these averages is well above the range of typical organic polymers (i.e., 1.3-1.5), which indicates a good choice of building blocks. As per the objectives of this work, we are able to identify a sizable number of compounds with RI values greater than 1.7, i.e., (a) 2851 (7.4%), (b) 22388 (13.1%), and (c) 9985 (99.8%). At or above the critical threshold of 1.8, we still find (a) 131 (0.3%), (b) 1252 (0.7%), and (c) 6961 (69.6%) compounds. The shift to considerably larger RI values also demonstrates the success of our strategy to prescreen the R 1 and R 2 individually and build PI candidates based on the top residues. (Details of the prescreening and screening results are provided in the Supplementary Material.) Fig. 5 shows the RI distribution of (a) R 1 , (b) R 2 , and (c) PI candidates containing each constituent building block (cf. the PI candidate library construction illustrated in Fig. 3). We find that R 1 and R 2 candidates containing building blocks B28 (anthracene), B25 (dibenzothiophene), and B24 (thianthrene) have the highest RI values, while those containing building blocks B17 (cyclopentadiene) and B16 (1,3,5-triazinane) show the lowest RI. The ranking of the building blocks is very similar for both R 1 and R 2 structures, which is unsurprising considering their general similarity. The PI candidates do not contain building blocks B10 (1,4-dithiane) and B14 (toluene with linking in 2,4 position), as these were missing in the top 100 structures of R 1 and R 2 that were used in their construction. The average RI values for all the other building blocks -with the exception of B26 (tetraphenylmethane) -are greater than 1.8. However, it should be noted that the number of PI candidates containing a specific building block is very variable. The red plot in Fig. 5 (c) shows the count of each building block in the 10,000 PI structures. The most common building blocks in the PI library are B1 (CH 2 -linker), B25, B2 (S-linker), B28, and B3 (O-linker), with B28 and B2 occurring in almost all PI candidates. Building blocks B1, B2, and B3 do not exhibit particularly large average RI values, however, given the construction template for the residues introduced in Fig. 3, they are statistically more likely to occur. The average RI value alone is thus not a sufficient metric to gauge the potential impact of each building block on the performance of PI candidates.</p><p>In the following, we analyze the contribution of the building blocks in the high-RI candidates of R 1 and R 2 , using the hypergeometric distribution analysis detailed in Sec. II D. We focus on the top 10% of the R 1 and R 2 candidates (with RI values greater than 1.687 and 1.711, respectively). We do not include PI candidates in this analysis as the selective construction scheme with its biased building block counts makes it less meaningful. Fig. 6 shows the resulting Z-scores, which identify the over-or underrepresentation of each building block in the high-RI R 1 and R 2 candidates compared to a random sample. The results point to B28 and B25 as the most prevalent building blocks in the top residues, and thus the most promising ones to consider for the design of high-RI polymers. This finding is in good agreement with the prior analysis of the RI averages and distribution. The ranking of the building blocks in both R 1 and R 2 are again largely the same, suggesting that the effect of the building blocks is similar for the somewhat different sequence in R 1 and R 2 .</p><p>The above Z-score analysis only yields insights into the pervasiveness of building blocks in the top 10% candidates. In order to gain a more comprehensive picture, we now evaluate the Z-score of each building block in each 10% segment of the R 1 and R 2 libraries. For this, we sort the R 1 and R 2 libraries by increasing RI value, divide them into 10 subsets, and perform a hypergeometric distribution analysis for each building block in each subset. The results are shown in Fig. 7. We observe certain general trends for the building blocks making up R 1 and/or R 2 :</p><p>1. The Z-score of some building blocks increases with increasing RI values, indicating a direct correlation (e.g., for B28, B25, B24, and B23).</p><p>2. The Z-score of some building blocks decreases with increasing RI values, indicating a negative correlation (e.g., for B17, B22, B10, and B16).</p><p>3. The Z-score of some building blocks goes through a maximum for intermediate RI values, indicating a corresponding correlation with average candidates (e.g., for B11, B26, B19, B18, B20, and B7).</p><p>4. The Z-score of some building blocks does not show a clear trend, indicating a lesser impact on and correlation with the RI values (e.g., for B1 and B6.)</p><p>In the bottom right corner of Fig. 7, we plot the average RI values of the R 1 and R 2 candidates in each of the 10 subsets. We observe that the R 2 structures have somewhat higher RI values in comparison to R 1 structures. The principal difference in R 1 and R 2 candidates is the number of aromatic building blocks, i.e., two vs three, and the higher content of these moieties correlates with higher RI values.</p><p>In addition to analysing the influence of individual building blocks on the RI values, we also study the potential impact of building block pairs. For this we calculate the joint Z-scores of all possible building block combinations in the top 10% candidates. The results for R 1 and R 2 are shown in Fig. 8. This analysis reveals some dependence on particular building block combinations. For instance, B23, B24, and B25 perform significantly worse when paired with B4 and B5, but exhibit large Zscores in combination with B2 and B3. However, overall we find the impact of individual building blocks to be the dominant factor. For instance, B28 has the largest positive Z-score, regardless of its counterpart.</p><p>Overall, we find anthracene (B28) to be the most promising moiety, a surprising finding and somewhat contradictory to our initial hypothesis considering that it does not contain sulfur. However, it is in line with (and independently confirms) other efforts in the community that aim to integrate anthracene into high-RI optical polymers (see e.g., Ref. [59]). The other outstanding moieties are the S-heterocycles dibenzothiophene (B25) and thianthrene (B24). Naphthalene (B23) shows some promise as well. The O-and S-linkers (B2 and B3, respectively) and to lesser degree SO 2 (B6) outperform the carbon-based linkers.</p><!><p>In this study, we demonstrated that the data-driven in silico approach that is being advanced by us and others can rapidly and efficiently assess the properties and performance potential of high-RI candidate compounds, identify numerous leads for next-generation PIs, and elucidate structure-property relationships that form the foundation for rational design rules. By combining our RI prediction model with virtual high-throughput screening techniques, we characterized candidates on a large scale at a fraction of the time and cost of traditional studies. We identified high-value building blocks (e.g., anthracene, dibenzothiophene, thianthrene) and structural patterns (e.g., S-and O-linkers) that correlate with large RI values. Correspondingly, we identified regions in chemical space in which we can hope to maximize the RI values of PIs. These guidelines allow us to target specific molecular motifs and create polymers with exceptional optical properties. In future experimen- tal work, we will utilize these guidelines and pursue the promising candidates that have emerged.</p><!><p>Electronic supplementary material accompanies this paper and is available through the journal website. It provides details of the computational data displayed in the figures throughout this paper or that were used in the statistical analysis. We also give detailed definitions of all statistical metrics used in this work.</p><!><p>The authors declare to have no competing financial interests.</p>