Systematic Computational Design and Optimization of Light Absorbing Dyes.

We present a workflow to aid the discovery of new dyes for the role of a photosensitive unit in the dye-sensitized photo-electrochemical cells (DS-PECs). New structures are generated in a fully automated way using the Compound Attachment Tool (CAT) introduced in this work. These structures are characterized with efficient approximate density functional theory (DFT) methods, and molecules with favorable optical properties are suggested for possible further use in DS-PECs. As around 2500 structures are generated in this work, and as we aim for still larger volumes of compounds to screen in subsequent applications, we have assessed the reliability of low-cost screening methods and show that simplified time-dependent density functional theory (sTDDFT) provides a satisfying accuracy/cost ratio. From the dyes considered, we propose a set that can be suitable for panchromatic sensitization of the photoelectrode in DS-PECs to further increase DS-PEC efficiency.

Introduction

Our development into a technologically advanced society has been accompanied by a rapid depletion of fossil fuels that makes the search for a sustainable energy solution of utmost importance. In the long term, one would like to harvest sunlight and use this energy to produce chemicals that can be used as fuels or feedstock. One way to achieve this is to split water inside a photo-electrochemical cell and obtain energy-rich molecular hydrogen.[1,2] Further processing could then lead to synthesis of methanol, ethanol, heavier hydrocarbons, and synthetic fuels.

Dye-sensitized photo-electrochemical cells (DS-PECs) are promising candidates for this efficient transformation of solar energy. However, the oxidative half-reaction that happens on the anode of the DS-PECs remains a challenging part yet to be fully optimized.[3,4] The interface between the electrode and photosensitive component is notoriously complex as it involves processes of light absorption, charge separation, and charge transport. It consists of a semiconductor sensitized with a light absorbing molecule (dye) that is coupled to a water oxidation catalyst (WOC), either via a covalent bond or coadsorbed on the semiconductor’s surface.[2,5] After absorbing a photon, the excited dye should promote charge separation and transfer an electron to the semiconductor. Separating these processes into different parts of the device, in analogy to natural photosynthesis,[6,7] allows individual optimization of the device units[8]—the dye and the semiconductor. In this work, we will focus on optimization of the dyes.

An optimal dye has characteristics such as a wide range of absorption for visible light, affinity toward fast electron injection into the semiconductor, and a low charge recombination rate. After electron injection, the oxidized dye must go back to the initial state by activating the WOC. In addition, the dye has to be chemically stable and has to show a strong anchoring to the semiconductor surface. The main advantages of employing metal-free organic dyes are their low cost and diversity.[9] There are many studies on organic dyes in many fields,[10] but we have not been able to find studies that systematically explore the large chemical space of these molecules. We believe that such studies, which are feasible with modern computational techniques, will enable us to quickly identify interesting new targets for synthesis and experimental validation of their properties.

In this paper, we present a tool for automated, systematic generation of molecular structures coupled to a workflow for computational screening of the generated molecules. The step into the area of big data requires automation of all the aspects of scientific research, from preparation to data generation and analysis. As we need to evaluate spectroscopic properties, the use of quantum chemical methods is mandatory, and the computational efficiency of the procedure is of importance. After generating structures by adding one or more substituents from a library to the core molecule, we therefore test the reliability of a number of approximate density functional theory (DFT) methods for prediction of the key absorption characteristics. We thereby focus on density functional tight-binding (DFTB),[11,12] time-dependent DFTB,[13] and simplified time-dependent density functional theory (sTDDFT),[14,15] as these methods are all capable of quickly evaluating the position and intensity of absorption peaks in the visible part of the solar spectrum. As a reference method for validation of their accuracy, we use time-dependent density functional theory (TDDFT) as well as comparison with experimentally known spectra.

We first apply our procedure to the 1,4,5,8-naphthalenediimide (NDI) dye as a significant number of its derivatives have been explored experimentally, providing good possibilities for validation of the theoretical predictions. After testing the procedure for this class of dye, we then apply the workflow to the perylene diimide (PDI) and perylene triimides (PTI)[16,17] as well, leading in total to almost 2500 derivatives that are evaluated.

Spectroscopic Properties of NDI

NDI is a small member of the family of aromatic diimides with a simple structure and relatively easy synthesis.[18,19] NDI is well-known for its application in artificial photosynthesis, and numerous scientific papers are proof of its extraordinary chemistry.[10,20−22] In this work, we focus on NDI’s light absorbing properties.

Non-substituted NDI absorbs in the UV region, and functionalizing of the naphthalene core (Figure ) can move the absorption to the visible region. Depending on the nature of ligands, electron withdrawing or electron donating, NDI’s absorption can be tuned throughout the solar spectrum, resulting in the so-called “rainbow” collection.[21] Substitution on the imide (Figure ) side has a negligible effect on the optical properties of NDI,[23] and it is often used to control solubility.

Figure 1

NDI molecule with highlighted possible substitution positions: core substitution (green), shoulder substitution (red), and substitution on the imide (blue).

Substitution with different types and numbers of functional groups has been addressed in the literature,[21,99,24] but most commonly, the effects of mono-substitution, 2,6-di-substitution, and 2,3,6,7-tetra-substitution (Figure ) are considered. It is shown that electron-donating functional groups like alkoxy,[25,26] alkylamino,[26] alkylthio,[25] and thiophenes[27] affect optical properties via electronic effects (mainly mesomeric and inductive) on the frontier molecular orbitals, causing a decrease of the HOMO/LUMO gap while increasing the HOMO and LUMO energies.[28,29] The lowest transition in NDI is susceptible to the substituent effect, while the higher transition is unaltered and resembles the properties of non-substituted NDI.[21,25]

As NDIs are planar π systems that tend to aggregate, an additional advantage of introducing the ligands is that the steric hindrance effect can be used to prevent the aggregation, making π–π stacking of the NDIs less favorable.

Methods

All computational methods used are included in the Amsterdam Modeling Suite (AMS) from Software for Chemistry and Materials (SCM).[30] Initial geometries for the core and the library of ligands are prepared with the SCM Graphical User Interface (GUI). Generation of new structures is done with the Compound Attachment Tool (CAT).[31] CAT is a Python code employing, among others, the Python Library for Automating Molecular Simulation (PLAMS)[32] for generating structures. The spectroscopic characterization is done with a workflow in which a structure optimization is followed by time-dependent (approximate) DFT calculations. The full workflow is depicted in Figure and can be divided into three parts: (i) generating new structures, (ii) calculating the properties, and (iii) processing and analyzing generated data. We will discuss these parts in more detail below.

Figure 2

Workflow scheme for the NDI example. The workflow is separated into three sections: (i) generating new structures—attaching the ligands to the substitution spots (highlighted in green) on the NDI in such a manner to avoid steric clashes; (ii) calculating the properties—running the geometry optimization and forwarding the optimized structure to the input for time-dependent calculations; and (iii) processing and analyzing generated data—manipulating data, visualizing data, and applying the selection criteria.

Structure Generation

The task of CAT is to make multi-substituted molecules by combining the core dye with a library of possible ligands. The user of this program needs to supply the sets of ligands and cores as a list of structure files in a suitable (xyz) format and needs to define the substitution spot(s) by selecting the hydrogen atom(s) to be substituted at the core and at the ligands. The coordinates of the hydrogen atom at the unsubstituted core serve as the initial position of the new substituent. Figure shows how molecule attachment works in CAT. Initially, we treat the ligand as a rigid body that is attached and rotated around the attachment bond axis to determine the position with minimal steric repulsion. After this rotation along the bond axis, the code will run an additional force field optimization in which the ligand is allowed to distort from its free molecule structure to minimize steric repulsion of the ligand with the rest of the molecule.

Figure 3

Visualization of CAT’s mono-substitution process; (1) rotation of the ligand in such a way that bonds of specified atoms are parallel; (2) translation of the ligand to the position of a specified atom at the core; (3) rotation around the bond axes in order to avoid steric clashes.

Multiple substitutions are implemented as series of individual substitutions on the specified substitution spots. As the dye molecules often possess symmetry, it is of importance to remove redundant structures by considering permutation symmetry. As we aim to have a procedure that will work for arbitrary input structures that may not always be generated by applying symmetry constraints in the structure optimization, we leave the specification of symmetry elements to the user. This also easily allows for the use of approximate symmetry, e.g., a larger symmetric dye that is functionalized in a way that breaks the exact symmetry. With this philosophy, the generation of unique derivatives may use techniques for permutational isomers.[33] CAT contains a function to eliminate the permutational isomers, taking into account the type of ligand, substitution sites on the core molecule, and symmetry type of the substitution site. For determining the number of unique derivatives, it thereby uses Polya’s enumeration theorem.[34,35] For the cases with up to four substituents that we consider in this work, this is primarily needed for the class of tetra-substituted molecules where the use of Polya’s theorem gives a significant reduction of the number of molecules that need to be considered explicitly.

A new single bond (between the core and ligand) increases the number of conformational isomers, in addition to conformational isomers of ligands and cores alone. The current version of CAT does not search for the lowest energy conformation. However, for a future version, we plan to include an interface to one of the already existing codes for conformer search.[36,37]

Prediction of Spectroscopic Properties

The initial molecular structures are optimized with the DFTB3[38] method using the 3-ob parameter set.[39−42] Convergence to a stable minimum is checked by calculating the lowest frequency normal modes.[43] We allow for some numerical noise and discard only structures with imaginary frequencies above 20 cm–1 in this step. We proceed with the remaining structures to the evaluation of molecular properties and compute the 10 lowest-lying singlet excited states with TDDFTB[13] with the same 3-ob parameter set, as well as sTDDFT[14] employing the range-separated[44] Coulomb-Attenuated B3LYP (CAM-B3LYP)[45] functional and the DZP basis set. The parameters for sTDDFT in combination with CAM-B3LYP were set as recommended in the paper by Risthaus, Hansen, and Grimme.[15] All the calculations are performed for isolated molecules in the gas phase.

For a subset of the final structures, properties are also computed with a reference TDDFT method, for which we also took the CAM-B3LYP[45] functional and the DZP basis set. CAM-B3LYP was reported to be reliable in predicting the low-lying excited states for the π-conjugated systems in general[46−48] and was recently favorably evaluated for NDIs in particular.[49]

Data Analysis

The results obtained from the calculations were processed using Python packages for data analysis and visualization, used in conjunction with SCM-GUI for visualization of molecular orbitals and production of simulated spectra.

Selection criteria in this analysis are the dye’s optical properties—excitation energies and oscillator strengths. Ideally, the lowest excited state should be above the threshold for water oxidation reaction that is about 950 nm (1.35 eV),[50] while the upper limit could be set at 3.2 eV due to the low intensity of UV light on the Earth’s surface. However, as we work with approximate DFT methods and apply additional approximations in our model (most importantly, the neglect of solvent and other environmental effects), we should also account for errors in the computed values when analyzing the results. For this reason, we discard molecules that have the lowest excitation with oscillator strengths below 0.001 as these are considered likely to be “spurious”.[51,52]

Of the remaining dyes, we consider the ones that have an intense, lowest transition above the energy threshold as the most promising. As not all promising molecules will be easy to synthesize, we also consider the use of a synthetic accessibility (SA)[53] score. This is a value that indicates how easy (1) or difficult (10) the synthesis will be.

Results

To generate the NDI derivatives, seven ligands with electron-donating characteristics were chosen from the literature.[21,25−27,54] The ligands contain sulfur, oxygen, or nitrogen and are numbered as follows: 1 is ethanethiol, 2 is thiophene, and 3 is bithiophene; 4 is ethanol; 5 is ethanamine, 6 is pyrrole, and 7 is 4-ethynyl-N,N-dimethylaniline.

The hydrogen atoms on the NDI core (Figure ) are substituted in the following order, the first position is 2 then 6, 3, and 7. If represented as a list, the first element of the list corresponds to substitution on position 2, the second element corresponds to substitution at position 6, etc.; if we use numbers to represent ligands as in Figure , we can represent NDI derivatives in a short notation as NDI(1,1) for 2,6-di-alkylthio-substituted NDI or even shorter as NDI-11. This way of substitution is found to be most common in the literature.[21,25,28,54] Using CAT, a total of 1015 NDI derivatives are generated: 7 mono-substituted, 28 di-substituted, 343 tri-substituted, and 637 tetra-substituted derivatives. The molecule sizes vary from 33 to 106 atoms. Experimental characterization is available for 11 molecules from this set. In order to assess the reliability of the predicted absorption spectra, we compared computed vertical excitation energies, EVE, with the absorption maxima, λmax, measured in solution (Figure ).

Figure 4

Ligands used to functionalize the NDI core. The atom highlighted in green is the connection point to the NDI core. 1, 2, and 3 are ligands containing sulfur (yellow); 4 contains oxygen (red); and 5, 6, and 7 include nitrogen atoms (blue). Carbon and hydrogen atoms are in gray and white, respectively.

Figure 5

Excitation energies of the isolated NDI derivatives calculated with TDDFTB, sTDDFT, and TDDFT compared to the experimentally obtained absorption maxima. Data points are sorted by the experimental values, increasing from left to right for each solution, chloroform (left) and dichloromethane (right).

The excitation energy values computed with the reference TDDFT(CAM-B3LYP) method show in general good agreement with trends in the experimentally obtained spectra. The method is able to predict the shifts observed in the experiment, e.g., the shift of the absorption band from two alkoxy substituents (NDI-44) to one alkoxy and one alkylamino substituent (NDI-54) and finally to two alkylamino substituents (NDI-55). Also, it predicts the correct electronic transitions; e.g., for the NDI-55, NDI-44, and NDI-54, the lowest energy peak is attributed to the HOMO→LUMO electronic transition, while the second peak is attributed to the HOMO-1→LUMO electronic transition. The second peak remains at the same position for these three molecules as it is a characteristic of the NDI core absorption.[25,26] An equivalent behavior occurs for tetra-substituted NDI-5555 compared to di-substituted NDI-55, with a shift toward shorter wavelengths.[25] On the other side of the energy scale, the NDI-4444 spectrum arises from three very close intense excitations, with the lowest on HOMO→LUMO transition followed by HOMO-1→LUMO and other higher transitions. The presence of the higher excitations has been suggested by Röger and Würthner,[25] as the mirror-image fluorescence band is absent. The molecule giving the lowest excitation energy in this figure (Figure ) is NDI-77, the NDI core di-substituted with the 4-ethynyl-N,N-dimethylaniline group. This large shift to lower energies is due to the extended π-conjugated system that is created by directly connecting the triple C–C bond to the NDI core and due to the electron-rich amino group at the other end of the ligand.[54]

Neglecting solvent effects on the geometry as well as on the spectra is partly responsible for the discrepancy between experimentally measured and calculated values. In the case of NDI-55 and NDI-5555, outliers in Figure , the automatically generated geometry differs from experimentally observed— the hydrogen attached to the nitrogen is not forming a hydrogen bond with the oxygen from the imide—and reduces the redshift.[26,28] For the non-polar chloroform solvent, we find an excellent quantitative prediction of the trend with TDDFT(CAM-B3LYP).

Turning now the attention to the performance of the approximate methods, we consider two possible screening methods, sTDDFT(CAM-B3LYP) and TDDFTB(3-ob), as well as the reference method, TDDFT(CAM-B3LYP). The statistical analysis is shown in Table for the same data points from Figure for values measured in dichloromethane. Looking at Table , we see the quantification of the consistent overestimation of excitation energies by TDDFT(CAM-B3LYP) and sTDDFT(CAM-B3LYP) that is already visible in Figure .

Table 1

Statistical Analysisa of the Computed Excitation Energies (eV) for the Isolated Molecule Compared with Experimental Results for the Molecule in Dichloromethane and Average Time (t) per Calculationb

method	MD	MAD	RMSD	R2	t
TDDFT(CAM-B3LYP)	0.58	0.58	0.59	0.92	27 min
sTDDFT(CAM-B3LYP)	0.39	0.41	0.41	0.88	4.5 min
TDDFTB(3-ob)	–0.17	0.25	0.30	0.60	3.5 s

MD stands for the mean deviation, MAD for the mean absolute deviation, RMSD is for the root-mean-square deviation, and R2 is the squared correlation.

Based on all molecules shown in Figure on two Intel Xeon nodes (48 cores) of the Dutch national supercomputer Cartesius.

The values for the mean deviation (MD) and mean absolute deviation (MAD) are higher than for TDDFTB(3-ob), which for some molecules underestimate and for others overestimate the excitation energy (MD = −0.17 and MAD = 0.25). However, the R2 value is more important as this quantifies whether shifts due to substitution are predicted correctly. Here, we see that, as expected, the correlation between theory and experiment is highest for the TDDFT(CAM-B3LYP), 0.92. sTDDFT(CAM-B3LYP) scores with 0.88, clearly better than TDDFTB(3-ob) that gives a value of only 0.60. The final column in Table is the average time per calculation, which is indicative of the computational cost of the workflow. Looking at its reliability and efficiency, sTDDFT(CAM-B3LYP) was found to be the optimal method for the purpose of massive screening of dyes. We also note that this method should be easier to combine with inclusion of solvent effects as it derives directly from a DFT ground state calculation. This analysis suggests that the range separation is indeed necessary, and instead of the local functional, TDDFTB(3-ob), its long-range corrected TDDFTB (lc-TDDFTB)[55] variant would be more appropriate for comparison. Unfortunately, this parameter set is not yet implemented in the SCM software suite.

Apart from solvent effects, the difference between computed values and experimental data can originate from various other factors.[56] A trivial empirical correction can be applied by shifting the excitation energy in dichloromethane down by the MAD value shown in Table . For sTDDFT(CAM-B3LYP), this would yield a modest 0.4 eV correction factor. To summarize the results of this first workflow application, we find that out of 1015 structures, 1013 are characterized as stable (no significant imaginary frequencies). Analysis of excitation energies gives 886 molecules with the lowest (non-spurious) excitations that fit the selection criteria. These excitations are shown in Figure and demonstrate that their energies span the energy range from 1.7 to 3.2 eV. By partitioning this range into smaller intervals and selecting the most intense dye from each interval (Table ), it is indeed possible to achieve panchromatic sensitization of the photoelectrode.

Figure 6

Absorption spectra based on the lowest excitation peaks from each promising compound predicted with sTDDFT. The color of the peaks resembles the color corresponding to the excitation. The absorption peak in black is from the non-substituted NDI (EVEsTDDFT = 3.56 eV; f = 0.39).

Table 2

Computed Lowest Excitation Energies with the Highest Oscillator Strength in Intervals of 0.2 eV between 1.7 and 3.2 eVa

energy range [eV]	NDI derivative	EVEsTDDFT[eV]	f	SA score
1.6–1.7	NDI-3333	1.68	0.53	7.91
1.7–1.8	NDI-3337	1.75	0.64	5.84
1.8–1.9	NDI-3377	1.83	0.82	6.87
1.9–2.0	NDI-3677	1.97	0.87	6.71
2.0–2.1	NDI-5677	2.05	1.02	7.95
2.1–2.2	NDI-5477	2.20	1.24	7.78
2.2–2.3	NDI-775	2.21	1.09	7.58
2.4–2.5	NDI-5447	2.43	0.66	7.38
2.5–2.6	NDI-126	2.53	0.46	7.13
2.6–2.8	NDI-5462	2.70	0.40	7.32
2.8–3.0	NDI-5454	2.83	0.35	7.05
3.0–3.2	NDI-444	3.10	0.27	6.80

Results are obtained with the sTDDFT method.

In Table , we notice that all compounds that have the highest absorption intensity contain ligand 7. This ligand gives rise to an extended π conjugation of the NDI core, which reduces the HOMO-LUMO gap of the NDI and thus lowers the excitation energy. The molecule with the highest oscillator strength has 4-ethynyl-N,N-dimethylaniline (7) at positions 2 and 6; with a third substituent, ethanamine (5), at position 3; the fourth substituent is ethanol (4) at position 7. Due to its low excitation energy and high oscillator strength, it can be considered as a promising molecule for further study. The derivatives NDI-3333 and NDI-77 have indeed already been recognized in the literature as promising candidates,[27,53] which confirms the workflow reliability. The last column in Table shows the SA score, which for all molecules falls between 6 and 8. This value is close to the value estimated for 11 experimentally known structures (Figure ) that are rated from 6.7 to 7.9. We also note that the SA score of non-substituted NDI is 6.44, which makes this value not very selective. This probably also has to do with the limited diversity of our ligands, and we believe that this value may become more useful for larger sets for which further reduction of the number of candidate structures is needed.

To proceed further into uncharted chemical space, we take three less well-known members of the perylene family and combine these with the same set of ligands. We hereby consider perylene diimide (PDI) and two core-extended PDIs, perylene triimides PTI1 and PTI2 (Figure ).

Figure 7

Three perylene-based chromophores (a) PDI, (b) PTI1, and (c) PTI2. Highlighted in green are substitution spots.

The PDI core can be substituted at 4 positions (Figure a), giving rise to the same number (1015) of mono-, di-, tri-, and tetra- substituted derivatives as the NDIs discussed previously. Calculated with sTDDFT, the lowest excitation energy of non-substituted PDI is 2.43 eV with an oscillator strength of 0.77. Compared to non-substituted NDI, its absorption is shifted dramatically toward the lower energies due to the vertical extension of the naphthalene core. Electron-donating substituents additionally shift the absorption to lower energies, but the effect of substitution is less diverse compared to the NDI core. PDI derivatives have the lowest excitations from 1.6 to 2.3 eV. Nine hundred eighty-six molecules had a stable structure from which 488 fit to the criteria (Figure ), as shown together with non-substituted PDI (black) in Figure . We see a systematic decrease in oscillator strength at lower energies and list in Table the most promising PDI derivatives for each energy range.

Figure 8

Absorption spectra based on the lowest excitation peaks from each promising PDI derivative calculated with sTDDFT. The color of the peaks resembles the color corresponding to the excitation. The absorption peak in black is from non-substituted PDI.

Table 3

Computed Lowest Excitation Energies with the Highest Oscillator Strength in Small eV Intervals in the Energy Range of PDI Derivativesa

energy range [eV]	PDI derivative	EVEsTDDFT[eV]	f	SA score
1.6–1.8	PDI-7133	1.77	0.39	6.96
1.8–1.9	PDI-5512	1.90	0.46	7.61
1.9–2.0	PDI-6744	1.93	0.57	7.92
2.0–2.1	PDI-4444	2.06	0.66	7.34
2.1–2.2	PDI-444	2.14	0.69	7.16
2.2–2.3	PDI-44	2.22	0.73	6.99
2.3–2.4	PDI-4	2.33	0.76	6.83

Results are obtained with sTDDFT(CAM-B3LYP)

The PTI1 core has three substitution spots, one at the core-extending imide group and two equivalent positions on other side of the core (Figure b). For this molecule, we first attach ligands to the extended side of the core leading to seven possible mono-substituted molecules (or new cores). Ligands are then attached to the two equivalent spots on the left side of the core, leading to 252 molecules in total, with 35 derivatives for each new core. The PTI2 core has only two positions for substitution (Figure c), which therefore leads to 35 PTI2 derivatives, seven mono-substituted and 28 di-substituted molecules. None of these 35 molecules have imaginary frequencies. The non-substituted PTI1 and PTI2 have almost the same lowest excitation energy, 2.88 and 2.89 eV, respectively, with also the same oscillator strength of 0.67. This indicates that the ligand attached on the core-extending imide group has no effect on the absorption properties.

For these molecules also, other ligands do not affect the optical properties significantly if attached to the core-extending imide group. In Figure , we noticed clusters of excitation peaks in the set of promising derivatives. These are four derivatives where X is a ligand (1, 4, 5, and 6). Y and Z are ligands that determine the position of the peak, and the core-extending imide group is acting as a node. Other three ligands at X are showing the same peaks, but the criterion that the lowest excitation is the most intense one is not fulfilled. Figure (bottom) shows HOMO and LUMO energies for the promising derivatives. The black squares show the HOMO and LUMO energies for the non-substituted PTI1. Compared to non-substituted PTI1, the LUMO energies are slightly higher, up to 0.5 eV. The effect on the HOMO level is a shift of up to 1.5 eV. The lowest energy peak at 2.3 eV (Figure ) belongs to PTI1-7-6. Table shows the most promising PTI1 derivatives for each energy range. PTI2 provides six promising structures (Table ), with the effects of its ligands being quite similar to PTI1.

Figure 9

(top) sTDDFT(CAM-B3LYP) prediction of absorption spectra based on the lowest excitation peaks from each promising PTI1. The color of the peaks resembles the color corresponding to the excitation. The absorption peak in black is from the non-substituted PTI1 core; (bottom) HOMO and LUMO energies in eV for the promising PTI1 derivatives.

Table 4

Lowest Excitation, Oscillator Strength, and SA Score for the Promising PTI1 Derivatives

energy range [eV]	PTI1 derivative	EVEsTDDFT[eV]	f	SA score
2.3–2.4	PTI1-66	2.38	0.44	7.39
2.4–2.5	PTI1-65	2.50	0.44	7.27
2.5–2.6	PTI1-12	2.54	0.40	7.24
2.6–2.7	PTI1-644	2.67	0.58	7.42
2.7–2.8	PTI1-644	2.72	0.59	7.23
2.8–2.9	PTI1-1	2.90	0.70	7.06

Table 5

Lowest Excitation, Oscillator Strength, and SA Score for All the Promising PTI2 Derivatives

PTI2 derivative	EVEsTDDFT[eV]	f	SA score
PTI2-6	2.36	0.44	7.37
PTI2-5	2.47	0.42	7.25
PTI2-2	2.51	0.38	7.36
PTI2-1	2.62	0.42	7.37
PTI2-44	2.64	0.54	7.40
PTI2-4	2.70	0.55	7.20

For the application in DS-PEC dyes, besides the strong light harvesting, the presence of a potential gradient is also important. For an example of a triad system, as in work of Monti, De Groot, and Buda,[57] the light absorbing dye also needs to fit the molecular redox properties, ground state oxidation potential (GSOP), and excited state oxidation potential (ESOP) of the two other components of the triad. Calculating these quantities is better done with a dedicated procedure,[58] the performance of which we plan to investigate in a separate paper. We note that, on the basis of the data calculated with the current workflow, one may give a simple estimate of the GSOP/ESOP using Koopman’s theorem as shown in Table . The MAD of the LUMO estimate is high, but results are considerably more consistent compared to the HOMO energies and possibly useful to provide insights into the relative ESOP of one derivative as compared to the others. The better performance of the LUMO estimate of ESOP can also be used to obtain a slightly better estimation for GSOP by subtracting the computed excitation energy from the LUMO, GSOPest = E(LUMO) – EVE.

Table 6

Statistical Analysisa of the Simple DFT Estimates with the Experimentally Obtained GSOP and ESOP

method	MD	MAD	RMSD	R2
ESOP(LUMO)	0.79	0.79	0.80	0.76
GSOP(HOMO)b	–1.95	1.95	1.97	0.57
GSOP(HOMO)est, c	0.15	0.24	0.28	0.76

MD stands for the mean deviation, MAD for the mean absolute deviation, RMSD is for the root-mean-square deviation, and R2 is the squared correlation.

HOMO energy from DFT; bHOMOest = LUMO – EVEsTDDFT.

Conclusions

We presented a tool for generating derivatives by systematic functionalization of a molecule to aid in exploring the chemical space of molecular derivatives. This generation tool can be combined with workflows for the screening of large numbers of molecules on the basis of computed molecular properties that fulfill certain criteria. We applied such a workflow to find molecules with potential to be used as photosensitizers in DS-PECs, using desirable optical properties as selection criteria. Results are shown for NDI’s, PDI’s, PTI1’s, and PTI2’s derivatives. The derivatives create a very diverse set of chromophores allowing panchromatic sensitization of the photoelectrode if adsorbed together on a semiconductor surface. The outcome is around 1400 derivatives that fulfill the criteria to be used as photosensitizers in DS-PECs. The workflow applied to NDI selects already known sensitizers out of 1015 structures, which shows its reliability. We also considered estimating the GSOP and ESOP by HOMO/LUMO energy levels but observed that this simple procedure based on DFT(CAM-B3LYP) values is not predictive enough for this purpose. The CAT code is available as an open source and is thereby extendable to other applications, where systematic screening of compound derivatives is needed, like in the discovery of new photoswitches or in the discovery of new water oxidation catalysts.