Quantitative Structure-Property Relationships for the Electronic Properties of Polycyclic Aromatic Hydrocarbons.

This study presents a development in quantitative structure-property relationships (QSPRs) for research in organic semiconductor materials by introducing a new structural descriptor called "degree of π-orbital overlap" based on two-dimensional structure information of aromatic molecules. Application of this method to predict the electronic properties of polycyclic aromatic hydrocarbon (PAH) molecules, which are known to be the core component of many organic semiconductor materials, is presented. Results demonstrated that QSPRs based on the new descriptor can predict rather accurate band gaps, ionization potentials and electron affinities for a large number of PAHs compared to those explicitly calculated by density functional theory method. This research opens new possibilities for developing QSPRs for other organic semiconductor classes in future.

Introduction

Organic semiconductor materials have opened a new era for developments of new organic-based electronic technologies and are predicted to be an alternative to amorphous silicon-based materials.[1−5] There are four main groups of organic semiconductor materials: organic light-emitting diodes (OLEDs), organic field effect transistors, organic thin-film transistors, and organic photovoltaics, each of which has different physical properties.

Organic semiconductor materials for fabricating electronic devices must be stable toward oxidation agents (e.g., oxygen) and have sufficiently high charge mobility as well as good photoelectronic properties.[5,6] To screen for materials satisfying such requirements, molecular properties, namely, lowest excitation energy, ionization potential (IP), and electron affinity (EA) are often used. These properties can be normally predicted from accurate first-principles quantum chemistry calculations, such as density functional theory (DFT). However, it is expensive in both time and computational resources. A number of more cost-effective methodologies based on quantum chemistry framework have been developed for predicting the electronic properties of organic semiconductor materials. For example, a previous study[7] suggested that semiempirical molecular orbital (MO) could be used to predict trends in band gaps for molecules within a given class since these properties differed from those calculated by more accurate density functional theory (DFT) by only shifting constants. Another approach is to couple a “solvation” model with DFT to include condensed-phase effects in predicting the electronic properties of organic solids.[8] A more common computational approach is the use of quantum mechanical (QM) and molecular mechanical models for studying the characters of disordered organic semiconductors.[9] Regarding applications, a high-throughput study to quantify the differences in band gaps of isolated polycyclic aromatic hydrocarbon (PAH) molecules and its crystals was reported.[10] Also, static polarizability and ionization potential were used as primary descriptors correlating with highest occupied molecular orbital (HOMO)–lowest unoccupied molecular orbital (LUMO) gaps of halogenobenzenes.[11]

Nevertheless, all of these methods are QM-based and do require time and efforts to compute properties of each compound and thus are not desirable at the initial designing stage where experimentalists only need reasonable estimates for a large number of possible candidates to arrive at a short list from which they can be further ranked by more accurate quantum chemistry calculations prior to carrying out actual experiments. Unfortunately, there is no possibility to do so.

Data analytics methodology, such as “quantitative structure–property relationship” (QSPR),[12] has been widely applied in drug design for screening biological activity, physicochemical properties, and toxicological responses. QSPR provides empirical relationships between chemical structural properties and the physical, chemical, or biological ones, which may require extensive calculations or experiments. There are some QSPRs currently available for predicting properties of organic semiconductor materials, such as for glass-transition temperature of OLED materials,[13] for charge-transfer rates of polycyclic aromatic hydrocarbons (PAHs),[14] or for predicting the toxicity of PAHs.[15] It is well known that PAH is the core component of many organic semiconductor materials in actual applications. For instance, PAH is the main component for synthesizing thienoacenes and tetraarylethenes in fabricating OLED materials.[3,16−18]

To the best of our knowledge, there have not been any QSPRs available for predicting important electronic properties, namely, band gaps, ionization potentials, and electron affinities, for a general PAH class of molecules. The objective of this study is to develop QSPRs for predicting the electronic properties of PAH class of molecules. Such QSPRs would greatly assist in the design of new PAH-based organic semiconductor materials.

Computational Details

QSPR Method

The QSPR method provides a framework for developing mathematical correlations between the chemical structural features and its corresponding properties. The essential aspect of the QSPR technique is the exploration of property features encoded within its structural descriptors. The basis formalism is response = f(descriptors-chemical structure, physicochemical property). Empirical mathematical relationships can be derived from a given set of data (training set), and then the relations can be used to predict/estimate properties of any similar molecules not in the training set. The main challenge of this method is to arrive at a set of descriptors for a given set of properties. A good QSPR is the one with the least number of descriptors while having a wide range of applicability, namely, a large molecular space. For that reason, descriptors for construction of reliable QSPR model should have the following features:[12]

be relevant to a broad class of compounds,

be fast in its calculation and independent of experimental/calculated properties, and

have different values for structurally dissimilar molecules, even if the structural differences are small.

Hence, constructing a good descriptor is the key to the success of the QSPR approach.

Various categories have been used to describe a QSPR descriptor, namely, constitutional: the relative numbers of various atom types; topological: describing bond connectivity of atoms in the molecule; physicochemical: relating to its physical or chemical properties, such as solubility in a given solvent (typically water or lipid), dipole moment, and formal charge; structural: describing the three-dimensional (3D) size, shape, and surface properties of the molecule; quantum chemical: relating to properties from quantum calculations, such as partial charges, polarizabilities, multipole moment, etc.[19] The choice of certain descriptors following the criteria suggested above is important for developing good QSPRs.

Previously, a number of indices were introduced to provide a quantitative understanding of aromaticity, such as harmonic oscillator model of aromaticity (HOMA), nuclear-independent chemical shift (NICS), para-delocalization index (PDI), aromatic fluctuation index (FLU), and multicenter index (MCI).[20] Although some interesting correlations between these indices and properties of PAH molecules were found, the use of these indices for QSPRs is limited. For example, HOMA presented excellent linear correlations with bond energy, but it is for a single ring not for a general PAH class while NICS index for a group of linear and angular PAH.[21] Furthermore, these indices are QM-based descriptors. For the purpose of this study, a structural or topological descriptor should be used.

BG(1) Descriptor

To understand aromaticity in aromatic molecules, Clar and co-workers suggested a method to quantify resonance structures known as aromatic sextet patterns. According to Clar’s theory, PAH isomer with more sextet patterns is more stable. Although the original theory only provides a qualitative understanding of aromaticity, it motivated a number of application studies on different classes of aromatic compounds as well as further theory developments.[22−34] One of such studies is the development of aromatic sextet polynomials by Hosoya et al.[27] that quantifies aromaticity of PAH structures.[22−24,26−28,31,33,34] In fact, Hosoya and Yamaguchi[35] introduced the “sextet polynomial” to quantify the level of resonance in an aromatic molecule by classifying different resonance structures in terms of the resonance sextet number r(G, k), which is the number of Kekule structures from the G molecule that have k sextets (isolated benzene rings with three double bonds). By definition, r(G, 0) is unity for all cases. Also, r(G, 1) is equal to the number of hexagons in the molecule.

For instance, dibenz(a,c)anthracene has five Kekule structures having two sextets (at hexagons (1, 4), (1, 5), (2, 4), (2, 5), (4, 5) in Figure ), i.e., r(G, 2) = 5, and two structures having three sextet patterns, r(G, 3) = 2.

Figure 1

Dibenz(a,c)anthracene.

The sextet polynomial BG(x) is then defined aswhere m is the largest number of sextets in the molecule.

Dibenz(a,c)anthracene’s sextet polynomial is 1 + 5x + 5x2 + 2x3. When x = 1, BG(1) is the sum of the coefficients of the polynomial, which is the same as the number of Kekule structures of the aromatic hydrocarbon.

In accordance with Clar, the Kekule structure with the maximum number of sextet pattern represents the electronic ground state of a PAH, and the more sextets a PAH has, the more stable it is. This suggests that there may be a correlation between the BG(1) and the relative stability of PAHs. Nonetheless, there has not been any study on the relationships between BG(1) and the physical or chemical properties of PAH molecules. Exploring the possibility of using BG(1) as a topological descriptor is an objective of this QSPR study.

Degree of π-Orbital Overlap Descriptor

Although BG(1) is based on well-tested Clar’s theory, evaluation of its value is not simple and thus its use for QSPR development and application would be limited. In this study, we propose a new structural descriptor called degree of π-orbital overlap (DPO). It is based on a known fact that aromaticity of PAH molecule relates to the delocalization of π-orbitals on the planar structural framework. Such delocalization affects molecular orbital energies and hence ionization potential, electron affinity, and band gap properties. Furthermore, previous studies have shown that π-orbital delocalization in PAH is affected by its structural shape.[5,36]

There are four basic building blocks or four different ways for adding a benzene ring to an existing PAH for constructing a larger one (see Figure ). These blocks suggest four separate parameters, namely, a–d to describe DPO for a given molecule.

Figure 2

Building blocks for constructing a given PAH molecule; (a–d) are simplified models showing how these building blocks are related to the DPO a–d parameters.

Simplified cartoons shown in Figure suggest that one can think of π-electrons in PAH molecules as particle-in-a two-dimensional (2D) box with a more complex shape. DPO value, in this case, would be an effective length of a simplified 2D box.

Rule for Determining the DPO Value of a PAH Structure

Step 1. The first step is to determine the reference segment by applying the following rules in sequential order till a distinct reference segment is found.

Rules for reference segment.

The segment with the largest number of fused rings in the molecule is the reference segment.

If all segments have the same size, then the segment with the largest number of parallel fused C–C bonds orthogonal to its direction is the reference segment.

If rules (a) and (b) do not yield a unique segment, then the segment with the least number of overlayers is the reference segment.

For example, for the PAH molecule shown in Figure , segment 1 has the largest number of fused rings (rule (a)) and thus is selected to be the reference segment with the line of reference as shown.

Figure 3

Segment 1 is the reference segment.

In Figure , the top-row molecule has three segments with the same size and thus rule (a) cannot be applied. Using rule (b), structure I shows that the reference segment is preferred over structure II by having a larger number of parallel fused bonds orthogonal to the reference direction. In the bottom-row molecule in Figure , both rules (a) and (b) do not yield a unique segment. In this case, rule (c) helps to select the reference segment, as shown in structure III instead of structure IV.

Figure 4

Segments with arrows in I and III are chosen to be the reference segments.

Step 2. Assigning DPO value for each fused bond.

After selecting the reference segment, DPO values are assigned to each fused bond and depend on four structure-dependent parameters as described below.

Beginning with the reference segment, starting from the left most, each fused bond has a value successively of 1, 1 – a, 1 – 2a, 1 – 3a,... where a is a parameter. DPO values for fused bonds on overlayer segments that are paralleled to the reference one, as shown in Figure , are scaled by factors of d, where k is the order of its layer above or below the reference, i.e., the segment lays kth row above/below the reference and d is also a parameter.

Figure 5

Illustration of assigning DPO topological parameters a–d to fused bonds.

For a segment that forms a 120° angle with the reference one (or its fused bonds form a 60° angle with those of reference segment), the DPO value starts at the value b and is scaled by a factor of d for each subsequent fused bond, as shown in Figures and 6. Similarly, a segment forms a 60° angle with the reference segment (its fused bonds from a 120° angle with those of the reference segment) and DPO values start at the value c and are also scaled by a factor d for each subsequent fused bond. b and c are parameters. Examples of assigning DPO for these types of fused bonds are shown below.

Figure 6

Examples of assigning DPO values for type 2 fused bonds.

The total DPO value for the molecule is the sum of all DPO values assigned for all fused bonds.

Determination of DPO Parameters

Since a–d parameters describe different topological aspects of PAH structures, their values can be determined by a relatively small number of PAH molecules with specific structural features. These molecules along with their properties are listed in Table S1 in the Supporting Information. For instance, the parameter a can be determined by minimizing the root-mean-square difference (RMSD) errors in its correlation with band gaps (HOMO–LUMO gap) for linear acenes from three to eight fused rings, as shown in Figure a. Once the parameter a is determined, it can be used to determine parameters b and c using the structures in Figure b,c, in the same manner, respectively. Finally, with a, b, and c already known, d is determined using the same procedure with structures in Figure d.

Figure 7

Molecules used to determine DPO structural parameters. (a)–(d) are subgroups of PAH molecules that were used to determine the DPO parameters a, b, c, and d, respectively.

The optimal values of these structural parameters are a = 0.05, b = −1/4, c = +1/3, and d = 1/3.

We note that the DPO value is unitless as BG(1).

Illustrations

Examples for assigning DPO value for each fused bond and their total values are given below for two different PAH molecules. The total DPO values of compounds A and B are 2.28 and 2.48, respectively (Figure ).

Figure 8

(A, B) Two examples of how to assign DPO parameters to fused bonds.

Data Sets

In this study, two separate data sets were used, namely, the QSPR parameter development set consisting of all PAH molecules having four to six fused benzene rings, and the test set consisting of all PAH ones having seven and eight fused benzene rings. Structures and properties of the development set are given in Table S2a,b in the Supporting Information. The QSPR parameter development set is used to develop the linear QSPR parameters to correlate each electronic property with the selected descriptor. The test set is used to assess the accuracy of the developed QSPR in predicting corresponding properties.

Computational Methods

For each molecule in both data sets, its 3D structure is fully optimized at the PM7 semiempirical MO theory using the MOPAC2016[37] package and at the B3LYP/6-31+G(d) DFT method using the Gaussian 09 program. B3LYP/6-31+G(d) method is known to predict rather accurately in calculating frontier orbital energies.[38−40]

The PM7 method was used to verify our previous finding on the relative differences in its calculated HOMO–LUMO gap with those DFT calculations for molecules in the same class. The HOMO–LUMO gap is used to estimate the band gap.[11] According to Koopman’s theory, the ionization potential is estimated by the negative value of the HOMO orbital energy, and similarly, the electron affinity is estimated by the negative value of the LUMO one.

Results and Discussion

Assessment of BG(1) Descriptor

Figure A–C shows good linear correlations (all of the correlation coefficients R2 are from 0.85 to 0.99) of BG(1) values with band gaps (A); with ionization potentials (B); and with electron affinities (C) of isomers with the same number of fused rings. In Figure , it is clearly observed that the structure having larger BG(1) value is more stable by having a larger band gap. This is consistent with the suggestion by Clar and Ruiz-Morales, a PAH molecule having the higher “maximum number of resonant sextets”[27] is more stable compared to its isomers with fewer aromatic sextets.[23,27,36] This result verifies for the first time the conclusion of Clar and Ruiz-Morales[23,36] about the possibility of existing correlations between the aromatic sextet numbers and the electronic properties of PAHs.

Figure 9

Linear correlations of BG(1) with HOMO–LUMO gaps (A); with ionization potentials (B), and with electron affinities (C) calculated using the DFT method; blue diamond dots indicate four fused rings, red squares indicate five fused rings, and green triangles indicate six fused rings.

However, these correlations are limited to only isomers having the same number of fused rings. For a larger molecular PAH space with a different number of fused rings, there is no correlation found. In addition, BG(1) values are difficult to calculate for complicated structures having more than seven fused rings using Hosoya’s[27] rule. Consequently, BG(1) descriptor is not appropriate for developing QSPR.

Assessment of DPO Descriptor

Linear correlations between DPO values and the three electronic properties are shown in Figure A–C.

Figure 10

Linear correlations between DPO and HOMO–LUMO gaps (A), ionization potentials (B), and electron affinities (C) of PAHs from four to six fused rings calculated using the DFT method.

This figure shows excellent linear correlations between DPOs and the electronic properties of PAHs with correlation coefficients (R2) all greater than 0.9. In addition, these linear relationships can cover a wide range of PAH molecules having different numbers of fused rings. Equations for these correlations of band gaps, IP, and EA properties with DPO are given in Table . Comparing to those in Figure , these results suggest that DPO is much better than BG(1) for using as a topological descriptor in developing QSPR. Moreover, calculating DPO value for any PAH is rather straightforward as described above.

Table 1

QSPR Equations for Band Gap, Ionization Potential, and Electron Affinity of PAH (All in Electronvolt)

electronic property	DFT equation	shifted PM7 equation	PM7 equation
band gap	y = −0.65x + 4.68	y = −0.51x + 4.44	y = −0.51x + 8.06
		c1 = 3.62
ionization potential	y = −0.30x + 6.04	y = −0.24x + 5.94	y = −0.24x + 8.75
		c2 = 2.81
electron affinity	y = 0.35x + 1.36	y = 0.26x + 1.50	y = 0.26x + 0.70
		c3 = −0.80

Scaling to DFT Accuracy

Semiempirical MO method, such as PM7, is much more cost-effective when calculating the electronic properties for a large number of molecules, compared to DFT methodology. Thus, it is more convenient to use PM7 to explore different QSPRs and to develop new descriptors with a condition that PM7 calculating properties can be scaled to those of DFT values. Our previous study found that PM6 band gap data could be shifted by a constant to match values calculated using the DFT method.[7] Here, we test if the same observation is applied for the PM7 method.

For each property, shifted PM7 values are calculated by adding the shifting constant c, which is defined as the average absolute difference between the PM7 and DFT values.

QSPRs for electronic properties from PM7, shifted PM7, and DFT methods are plotted in Figure , and their corresponding equations are listed in Table . The present results further confirmed our earlier finding that semiempirical MO method can correctly predict trends in the electronic properties, namely, band gap, IP, and EA, as its linear equations differ from those of DFT data mostly by shifted constants. Root-mean-square difference (RMSD) values between DFT values and shifted PM7 ones for band gap, ionization potential, and electron affinity, respectively, are 0.123, 0.052, and 0.076 eV. In addition, available experimental values[41,42] (triangles in Figure ) for the band gaps of some PAH structures are also plotted. The results illustrate that both shifted PM7 and DFT QSPRs are within the experimental uncertainty. This suggests that one can use the PM7 method to develop QSPRs or for screening for trends in the electronic properties of materials. Corrections can be made by calculating the shifting constants using a small number of molecules in the training set.

Figure 11

Comparisons between PM7, shifted PM7, and DFT results for (A) HOMO–LUMO gaps, (B) ionization potentials, and (C) electron affinities. The red dots indicate DFT calculations; the blue squares indicate PM7 calculations; the triangles indicate experimental data; the solid lines indicate DFT calculations; the dot lines indicate shifted PM7; and the dash lines indicate PM7 calculations.

All of the equations are shown in Table .

Assessment of DPO-Based QSPRs

To assess the accuracy and predictability of the DPO-based QSPRs for the electronic properties of PAH molecules, properties of all molecules in the test set were also calculated at the same DFT level theory and plotted against those predicted by QSPRs, as shown in Figures –14.

Figure 12

Plots of QSPR vs DFT for HOMO–LUMO gaps of all molecules in the test set (A) for shifted PM7 QSPR equation and (B) for DFT QSPR equation.

Figure 14

Plots of QSPR vs DFT for electron affinities of all molecules in the test set (A) for shifted PM7 QSPR equation and (B) for DFT QSPR equation.

Plots of QSPR vs DFT for ionization potentials of all molecules in the test set (A) for shifted PM7 QSPR equation and (B) for DFT QSPR equation.

According to the results from Figures –14, QSPRs from shifted PM7 and DFT predict equally well the band gap and ionization potential. Statistically, those from shifted PM7 QSPR perform slightly better. The R2 values of correlations range from 0.96 to 0.98, which demonstrates that the DPO is an excellent descriptor for developing QSPRs for the electronic properties of aromatic systems.[19]

How To Use DPO-Based QSPRs

For a given PAH molecule, the first step is to determine its DPO value from the rule described above: (a) choose the reference segment; (b) assign each fused C–C bond its value; and (c) sum all of these values to give DPO.

For band gap, use the QSPR y = −0.65x + 4.68 (where DPO is x); similarly, y = −0.30x + 6.04 for ionization potential equation, and y = 0.35x + 1.36 for electron affinity.

We note that higher DPO value predicts lower band gap for the PAH structure, which also means that the molecule is less aromatic.[11]

Issues with Planarity

We recall that the degree of π-orbital overlap descriptor presented in this study is based on the 2D topological information of the molecules. For PAHs, 2D structures/drawings are planar. Previous studies have reported a variety of structural characters that can affect the electronic properties of PAH, such as the maximum number of aromatic sextets and the holes in each Kekule resonance structure,[22−24,26−28,32−34,43] the zigzag or linear shape of the molecule structure,[36,44] as well as its structural planarity.[36,45] DPO can capture most except the planarity structural character. For that reason, one can expect nonplanar PAH molecules to be exceptions to the present DPO-based QSPRs. Below we provide several examples to illustrate this case.

According to DPO-based QSPRs, the band gaps of molecules shown in Table would decrease from left to right (1–3) from 2.65 to 2.36 eV due to more parallel fused bonds added to molecules, as suggested by the rule. Explicit DFT calculations for HOMO–LUMO, however, predict no significant difference between them. This is due to the fact that molecules 1–3 are nonplanar. Such nonplanarity, as shown in Figure , is caused by steric effects between the two segments. Due to such nonplanarity, the delocalization of electrons in π orbitals from one segment to the other is disrupted and thus band gaps of these molecules appear to depend only on the band gap of the longest segment. In principle, the degree of nonplanarity can be added as an additional correction descriptor to QSPRs although it would increase the complexity of the method and can be considered in a future study.

Table 2

DFT-Calculated and DPO-Based QSPR-Predicted Band Gap (eV) for Selected PAHs

Figure 15

Nonplanarity in molecule 3 in Table due to steric effects.

Conclusions

This study presents the development of new quantitative structure–property relationships to enable predictions of electronic properties, namely, band gaps, ionization potential, and electron affinities of PAH molecules without having to perform expensive explicit quantum chemistry calculations. We found that the structural property based on the well-known Clar’s theory of aromaticity provides good linear correlation with the electronic properties of PAH molecules, yet limited to only isomers having the same number of fused benzene rings.

A new descriptor called “degree of π-orbital overlap” based on 2D topological information of PAH molecules is developed. Results show that DPO-based QSPRs can accurately predict band gap, ionization potential, and electron affinity compared to those from explicit DFT calculations for a large class of PAH molecules. Comparing to previously reported aromatic indices, such as HOMA, PDI, FLU, MCI, and NICS, DPO establishes a much simpler and yet more general framework for correlating the electronic properties of PAH molecules with its 2D topological information. Consequently, it will open new possibilities for using DPO for further developments of QSPRs for organic semiconductor materials.