Jelena
Belić
^{a},
Arno
Förster
^{a},
Jan Paul
Menzel
^{b},
Francesco
Buda
^{b} and
Lucas
Visscher
*^{a}
^{a}Department of Chemistry and Pharmaceutical Sciences, Vrije Universiteit Amsterdam, De Boelelaan 1083, 1081 HV, Amsterdam, The Netherlands. E-mail: l.visscher@vu.nl
^{b}Leiden Institute of Chemistry, Leiden University, Einsteinweg 55, P.O. Box 9502, 2300 RA, Leiden, The Netherlands
First published on 1st December 2021
Sustainable solutions for hydrogen production, such as dye-sensitized photoelectrochemical cells (DS-PEC), rely on the fundamental properties of its components whose modularity allows for their separate investigation. In this work, we design and execute a high-throughput scheme to tune the ground state oxidation potential (GSOP) of perylene-type dyes by functionalizing them with different ligands. This allows us to identify promising candidates which can then be used to improve the cell's efficiency. First, we investigate the accuracy of different theoretical approaches by benchmarking them against experimentally determined GSOPs. We test different methods to calculate the vertical oxidation potential, including GW with different levels of self-consistency, Kohn–Sham (KS) orbital energies and total energy differences. We find that there is little difference in the performance of these methods. However, we show that it is crucial to take into account solvent effects as well as the structural relaxation of the dye after oxidation. Other thermodynamic contributions are negligible. Based on this benchmark, we decide on an optimal strategy, balancing computational cost and accuracy, to screen more than 1000 dyes and identify promising candidates which could be used to construct more robust DS-PECs.
In the present work we will focus on the dye on the photoanode of DS-PECs, for which the absorbed sunlight creates a photo-excited electron that is injected into the photo anode (usually TiO_{2}). The thus oxidized dye is restored to its initial state by accepting an electron from a water oxidation catalysts (WOC), which then starts the water oxidation cycle.^{7} To make this process efficient the dye has to fulfill several requirements. Prerequisite for the chain of electron transfers (from the dye to the semiconductor and from the WOC to the dye) to proceed is a suitable alignment of the dye's redox potentials with the neighbouring components: (1) the ground state oxidation potential (GSOP) has to be higher than the highest oxidation potential (HOP) of the WOC and (2) the excited state oxidation potential (ESOP) has to be lower than the edge of the anode's conduction band (CB) (Fig. 1). While these can be considered as minimal requirements, one may further narrow the search for optimal dyes by estimating that for reducing the energy loss and for an optimal rate of the electron transfer the potential of the ESOP state needs to be just slightly higher (∼0.3 eV)^{8} than the semiconductor's CB.
A further boundary condition for sustainable deployment of DS-PECs is the use of dyes made entirely from abundant material, such as fully organic dyes. Among these, perylene diimides (PDI) have been recognised as very promising for use in photovoltaic devices^{7,9} due to their favourable optical and electrochemical properties, as well as their stability, relatively cheap synthesis and tunability. They have found a place in different types of photo-electrochemical devices, both at the anode^{10–13} and at the cathode.^{14} Also other dyes related to the perylene diimide have proven to be efficient components of such systems, like naphtalene diimide (NDI)^{15} or perylene triimides (PTI).^{16}
Given this versatility of the perylene group of dyes, there are thousands of dyes that could be synthesized and tested for their performance. To better narrow down the search for optimal performers, the use of computational modelling is of interest. With a clear idea and well-defined objectives, such a computational screening of the chemical space can strongly reduce the amount of trial experiments. Important is thereby to take into account the uncertainties that the construction of model and the choice of theoretical methods bring. In a complex system there are many factors that affect the final performance and theory should merely narrow down the search rather than finding the perfect candidate. A great advantage is the modularity of the DS-PEC system, thus one may study its molecular components separately and define suitable components before constructing an overall model. Studies on molecular components therefore encompass dyes with favourable light absorption properties,^{17} anchoring groups between the dye and the semiconductor for ultrafast electron injection^{18} and the WOC–dye complex for driving the water splitting process.^{19,20} Analysis and screening of the individual components will help to define optimal DS-PECs and to improve upon the hydrogen production efficiency. In this work we focus on finding dyes that fit the system's thermodynamic requirements, in particular those that have a suitable GSOP.
The GSOP can be predicted by calculating the Gibbs free energy of the oxidation reaction,^{21} the difference between the solution-phase Gibbs free energies of a molecule in neutral and oxidized form. We can access these thermodynamic properties using generalized^{22,23} Kohn–Sham (KS)^{24} density functional theory (DFT).^{25} As reactions take place in a medium, solvent effects need to be accounted for which is most efficiently done with continuum solvation (CS) models.^{26} Furthermore we may distinguish between two major pathways that are commonly used to calculate solution-phase Gibbs free energy: use of a thermodynamic cycle (TC) and direct computation (DC) within a CS model.
Simpler than these two thermodynamical approaches is to approximate the GSOP by the vertical ionization energy of the solvated dye. Such a calculation can be done in three ways, one may compute the electronic energy difference by two separate KS-DFT calculations of the molecule in its two forms (before and after oxidation), or one may take the GSOP as the negative of the highest occupied molecular orbital (HOMO) energy of the neutral molecule. For the latter approach, it is important to take into account that semi-local or hybrid approximations to the exact functional of DFT do not describe the physics of addition or removal of electrons fully correctly.^{27–30} A theoretically more rigorous alternative is many-body perturbation theory (MBPT) in the GW approximation (GWA) to the electronic self-energy,^{31} where G denotes the single-particle Green's function and W is the screened Coulomb interaction. As in approaches based on total energy differences, the GWA accounts for the systems response to the hole created by ionization.^{32} It yields QP energies which can be directly identified with vertical electron addition and removal energies^{33} and is therefore increasingly used for energy level alignment in photo-catalytic interfaces^{34,35} or dye-sensitized solar cells.^{36–39}
The most commonly used way to validate theoretical predictions of the GSOP is a comparison to the oxidation potential measured by cyclic voltammetry (CV). The half-reaction (Dye → Dye^{+} + e_{anode}^{−}) taking place at the surface of the anode, with the dye redox couple [Dye/Dye^{+}] rapidly exchanging electrons with the electrode, is thereby considered to be electrochemically reversible. The measured half-wave potential E_{1/2} of the reaction is, at standard conditions, a very good estimate for the standard electrode potential: E_{1/2} ≈ E°,^{40} with E° being the driving force for electrochemical work.
As we are primarily interested in the GSOP here, the measured E_{1/2} potentials can directly be used as experimental reference data to validate the accuracies of the different methods to calculate the GSOPs. Treatment of the ESOP is a bit more involved as this depends on the GSOP as well as on the energy difference between the ground and excited state. The latter can be computed as the adiabatic energy difference between the ground state and the geometry-relaxed excited state, but also as a vertical excitation energy under the assumption that the electron transfer is much faster than the nuclear relaxation (Franck–Condon principle). Which picture is the most appropriate depends on the dye-semiconductor system; it has been observed for some semiconductors that the injection happens only after the relaxation.^{8} One may also note that for certain dyes the excited-state reorganisation energy is very small, making the differences between the schemes less important.^{41}
In this work we aim to construct and validate a fast procedure for the calculation of GSOPs which can be used for the automated screening of the GSOPs for different sets of dyes and their derivatives. In this construction step we consider accuracy as well as computational efficiency and discuss the effect of different approximations on these factors. This work represents a step towards full automation for the screening of dyes for use in DS-PECs – an automated strategy to obtain accurate energetics of the system's interfaces.
G^{i}_{sol}(g^{i}_{sol}) = E^{i}_{sol}(g^{i}_{sol}) + G^{i}_{therm}(g^{i}_{sol},T) | (1) |
The thermal contribution is the most demanding term in the calculation as it requires determination of the vibrational frequencies of the molecule. The G^{i}_{therm}(g^{i}_{sol},T) is computed within an ideal gas model and includes, besides zero-point vibrations, the translational and rotational contributions to nuclear kinetic energy. In the computation of the entropy term the quasi-rigid-rotor-harmonic-oscillator approximation model^{44} for weak (less than 20 cm^{−1}) vibrational modes is employed to avoid artifacts due to inaccuracies in these modes.
Given these individual Gibbs free energies, the GSOP can then be calculated as the difference (ΔG^{ox}_{sol}) between the solution-phase Gibbs free energies of a neutral and charged species (Fig. 2, bottom reaction, blue arrow). We will refer to this way of calculation as direct computation (DC). This procedure is rather straightforward and the final value of the adiabatic GSOP with DC and COSMO^{45} as CS model, ΔG^{DC}_{COSMO} is calculated using eqn (2):
ΔG^{DC}_{COSMO} = G_{sol}^{+}(g_{sol}^{+}) + G_{gas}(e^{−}) − G^{0}_{sol}(g^{0}_{sol}) | (2) |
An alternative is the use of a thermodynamic cycle (TC) method (Fig. 2, top reactions, red arrows). In the TC method, a gas-phase Gibbs free energy is augmented by solvation contributions to yield a solution-phase Gibbs free energy. For each oxidation state (neutral, 0, or ionized, +) three separate contributions are computed: the gas-phase Gibbs free energy G^{i}_{gas}(g^{i}_{gas}) (following eqn (1) for gas phase), the Gibbs free energy of solvation computed at the gas phase geometry g^{i}_{gas}, and the energy difference between the solvated and gas phase structure:
G^{i}_{sol}(g^{i}_{gas}) = G^{i}_{gas}(g^{i}_{gas}) + ΔG^{i}_{solv}(g^{i}_{gas}) + [E_{gas}(g^{i}_{sol}) − E_{gas}(g^{i}_{gas})] | (3) |
The equation for the adiabatic GSOP computed with the TC method and COSMO, ΔG^{TC}_{COSMO} is then similar to eqn (2) and reads:
ΔG^{TC}_{COSMO} = G_{sol}^{+}(g_{gas}^{+}) + G_{gas}(e^{−}) − G^{0}_{sol}(g^{0}_{gas}) | (4) |
The DC and TC strategies can be expected to give close agreement if the structure and vibrational frequencies of the dye are only weakly affected by solvation. If the stationary points in solution strongly differ from the ones in the gas phase, the DC and TC approach will deviate from each other.^{46,47} For the thermal correction the use of a thermodynamic cycle is then often considered to be more robust^{48} as CS models are parametrised to yield accurate Gibbs free energies of solvation (Fig. 2, vertical red lines).^{46,49–52} This implies that the thermal contributions calculated in a dielectric continuum in the (DC strategy) can contain corrections that are already incorporated in the CS model (via parametrization) which might result in double counting for some solvent effects. Ho^{47} compared the thermal contributions (particularly vibrational corrections) obtained by the two strategies, TC and DC, and showed that mean absolute errors in the redox potential for the two strategies differ by 10 mV from each other. In this work he also considered introducing a higher-level electronic structure theory, but found that this did not necessarily lead to better agreement with experiment, indicating that the description of solvation via a CS model is the main source of errors. Regarding the choice between the DC and TC approach, we furthermore note that in the TC approach evaluation of multiple solvents is more economical as the most demanding step of the calculation, determination of vibrational frequencies, is independent of the chosen solvent. If only one solvent is to be considered, the TC and DC strategies are similar in terms of computational costs. In DC we avoid the separate calculations for the gas and solution phase that are needed in the TC approach, but the effort needed for frequency calculations in solution is typically higher than in gas phase.
As a cheaper alternative in the TC pathway we can also employ the conductor-like screening model for realistic solvents (COSMO-RS),^{53,54} which performs statistical thermodynamics for the molecular surface interactions after the quantum chemical calculation. While COSMO provides the Gibbs free energy of solvation only (vertical red reaction), COSMO-RS provides the chemical potential of a solvated molecule and allows calculating all fluid phase equilibrium thermodynamic properties which includes the solution-phase Gibbs free energy as well.^{50,55} The low cost of this model lies in statistical thermodynamics part whose cost is negligible compared to quantum chemical calculations and which makes the explicit calculation of thermal contributions redundant, as it directly provides the solution-phase Gibbs free energy. The solution-phase Gibbs free energy thus reduces to a sum of the gas-phase energy and a COSMO-RS Gibbs free energy correction. The geometry changes due to solvation are ignored,
G^{i}_{sol,CRS}(g^{i}_{gas}) = E^{i}_{gas}(g^{i}_{gas}) + ΔG^{i}_{CRS,solv}(g^{i}_{gas}) | (5) |
Finally, to be able to compute a large number of molecules, we can reduce the time spent for geometry optimization. A workflow combining semiempirical quantum mechanical (SQM) techniques for geometry optimisations and frequency calculations and DFT electronic energy with COSMO is shown to provide reliable predictions of redox potentials.^{56,57} Similarly, we can employ the SQM method instead of DFT for geometry optimisation and combine it with DFT electronic energies and the COSMO-RS solvation model. The combination of SQM methods with a COSMO-RS thermodynamic calculation is very efficient and the limiting step in this workflow will be the single point DFT calculation including COSMO. A dedicated COSMO calculation is needed to define the molecule's reference state that is used in further COSMO-RS statistical thermodynamics calculations. As COSMO-RS provides the solution-phase Gibbs free energy, no frequency calculations are needed. We will denote the result of this composite method as ΔG^{screening}_{COSMO-RS}.
ΔE^{ox} = E_{sol}^{+}(g^{0}_{sol}) − E^{0}_{sol}(g^{0}_{sol}) | (6) |
(7) |
W_{0}(ω) = [ε_{RPA}(ω)]^{−1}v_{c}, ε_{RPA} = 1 + v_{c}π_{0}(ω), | (8) |
Σ_{ii}(ω_{i}) = ω_{i} − ε_{i}, | (9) |
Each substituent is given a single character for identification to shorten the notation for the complete dye in which the substituents are attached to the aromatic core, as shown in Fig. 4. The substituents containing electron donating atoms and groups are numbered from 1 to 7. The substituents 8, 9 and 0 are introduced only to enlarge the experimental data set available for validation and they are not used in the final screening procedure. As shown in Fig. 4, side chains of the substituents and cores were sometimes shortened for computational convenience. We thereby checked that the consequence of such replacements on the computed GSOPs are minor, and moreover, since this is done systematically for all the dyes this should give at most a consistent shift in oxidation potential value (ESI,† S1). For comparison of electrochemical properties from different experiments or for comparison of computational results with experiment we need to convert all values to a common scale, the absolute scale (against vacuum). Choosing Fermi–Dirac statistics as statistical mechanical formalism for the electron we set the absolute potentials for F_{c}/F_{c}^{+} couple as 4.98 V and for SCE as 4.52 V (ESI,† S1).^{43} The complete set of dyes with experimental oxidation potentials can be found in the ESI,† (Table S1).
To quantify the predicting ability of the methods we will use mean deviations and squared correlation – R^{2}. The magnitude of mean deviations mainly depends on the choice of the formalism to relate the reference electrode scale to the absolute scale. This choice might yield higher or lower mean deviations, but the shift will be consistent throughout the data set. For this reason we will rely on the R^{2} – strength of the linear correlation – for the evaluation of the performance. Having in mind possible errors from experimental and theoretical methods we monitored how sensitive R^{2} is to small variations in the reference values. Considering the errors that might occur as a consequence of side chain adjustments we tested how much R^{2} deviates from its initial value if we introduce a random error in the (+0.05, −0.05) eV interval. We found that the mean absolute deviation of the R^{2} converges to 0.01 (ESI,† S2). In conclusion, the methods for which the correlation with the experiment differs in ±0.01 would be considered to perform equally well.
The quality of density fitting and numerical integration were set to Good, except for the part with Hartree Fock exchange where quality was set to Normal. The SQM technique for GO in the screening workflow was GFN1-xTB.^{89} For the dyes that showed too large imaginary frequencies the GO calculation was repeated after tightening the convergence criteria for energy and nuclear gradients by a factor of 10 relative to the default values. The frequency calculations for hybrids were done numerically; their outputs contain all the thermodynamic properties, electronic bonding energy and nuclear kinetic energies, at room temperature. If present, negative frequencies were re-scanned to assure that they were not spurious (tolerating frequencies above −20 cm^{−1}).
Solvation effects were incorporated with the COSMO^{90} or COSMO-RS^{91} models, at a pressure of 1 atm and temperature of 298.15 K. All the calculations including the COSMO model were performed in dichloromethane, with a dielectric constant of 8.9. The default Amsterdam density functional (ADF) 2019.302^{92} settings were used: the cavity construction is Delley type,^{93} atomic radii are the corresponding van der Waals radii from the MM3 method by Allinger^{94} (with widely accepted increase of 120%^{54}). COSMO-RS requires an input from a COSMO calculation that defines the reference state which is the molecule in a perfect conductor, with infinite dielectric constant. In the reference COSMO calculation, the cavity construction is of Delley type and atomic radii are the Klamt atomic radii.^{54} The subsequent statistical thermodynamic calculation describes the effect of dichloromethane as a solvent.
The GW calculations were based on the solvated geometries optimized with B3LYP/TZ2P. We considered different variants of non self-consistent GW (G_{0}W_{0}) and partially self-consistent GW: G_{0}W_{0}@PBE0(HF = 0.4) (G_{0}W_{0} performed on top of a PBE0^{95,96} calculations with 40% exact exchange), G_{0}W_{0}@LRC-ωPBEH,^{97} eigenvalue-only self-consistent GW (evGW) based on the same two functionals (evGW@PBE0(HF = 0.4) and evGW@LRC-ωPBEH), QP self-consistent GW (qs GW) as well as QP self-consistent GW with self-consistency only in G, but not in W based on a PBE starting point (qsGW_{0}@PBE). All GW calculations were performed using the TZ2P and TZ3P^{71} basis sets. For G_{0}W_{0}@PBE(HF = 0.4) we additionally used the QZ6P basis set^{71} and extrapolated the IP to the CBS limit using
(10) |
The GW implementation in ADF follows the space-time method and is outlined in ref. 67 and 71 as well as in ref. 68 for qsGW_{0} and qsGW and the numerical settings chosen for all the GW calculations in this work follow the recommendations therein: We used 20 grid points in imaginary time and imaginary frequency each. For the G_{0}W_{0} and evGW we set NumericalQuality Good and for the qsGW_{0} calculations we used veryGood quality combined with the Good fit sets and thresholds for the HF part of the gKS and for the GW part of the calculations. For reasons outlined in ref. 67 we set Dependency Bas = 5 × 10^{−4} in the AMS input for G_{0}W_{0} and evGW calculations and Dependency Bas = 5 × 10^{−3} for qsGW_{0}.
Approach | Calculations^{a} | Method | Elapsed (step)^{b} | Total (core h)^{b} |
---|---|---|---|---|
a GO denotes geometry optimization of the neutral (n) or oxidized (o) molecule in vacuum (v) or in solution (s); f denotes frequencies; SP denotes single point calculation using the preceding geometry. All geometry calculations are started from the same initial geometry. b All calculations have been performed on a single 2.2 GHz intel Xeon (E5-2650 v4) node (broadwell architecture) with 24 cores and 128 GB RAM. c The TZ2P basis set and B3LYP-D3(BJ) functional has been used for all DFT calculations. Good numerical quality has been used throughout except for the Hartree–Fock part, where the Normal fit set has been used. d (QZ6P + TZ3P). | ||||
ΔG^{DC}_{COSMO} | GO(n,s) + f | DFT | 390 | 1066 |
GO(o,s) + f | 676 | |||
ΔG^{TC}_{COSMO} | GO(n,v) + f | DFT | 381 | 1109 |
+SP(n,s) | 1.4 | |||
GO(n,s) | 9.6 | |||
+SP(n,v) | 1.3 | |||
GO(o,v) + f | 692 | |||
+SP(o,s) | 2.5 | |||
+SP(n,v) | 1.3 | |||
GO(o,s) | 18.3 | |||
+SP(n,v) | 1.3 | |||
ΔG^{TC}_{COSMO-RS} | GO(n,v) | DFT | 9.6 | 32.7 |
+SP(n,s) | 1.4 | |||
GO(o,v) | 19.2 | |||
+SP(o,s) | 2.5 | |||
ΔE^{DFT}_{ox} | GO(n,s) | DFT | 9.6 | 12.1 |
+SP(o,s) | 2.5 | |||
ε ^{DFT}_{HOMO} | GO(n,s) | DFT | 9.6 | 9.6 |
ε ^{ GW }_{HOMO} | GO(n,s) | DFT | 9.6 | 47.7 |
+SP(o,s) | DFT | 2.5 | ||
+SP(n,v) | G _{0} W _{0} ^{ } | 29 + 6.6 | ||
ΔG^{screening}_{COSMO-RS} | GO(n,v) | GFN-xTB | 0.1 | 4.1 |
+SP(n,s) | DFT | 1.4 | ||
GO(o,v) | GFN-xTB | 0.1 | ||
+SP(o,s) | DFT | 2.5 |
As the GSOP represents the difference between the two different oxidation states of the same molecule, some contributions will cancel out. Here we elaborate on the effective contributions to the GSOP originating from solvation, geometry relaxation due to oxidation and due to solvation and thermal contribution. The energy contribution due to geometry relaxation related to the solvation in the TC method is added separately from electronic solvation effects, as shown in eqn (3).
We found that this energy difference for both, neutral or oxidized molecules is negligible, being on average 0.02 eV and independent from type or number of ligands (ESI,† S5). Therefore, for this type of molecules solvated in dichloromethane the energy contribution coming from the change in geometry due to the solvation is negligible. In the TC method, this means that evaluation of these contributions can be omitted to reduce the total computational cost of the ΔG^{TC}_{COSMO} strategy.
On the other hand, the electronic solvent contributions are crucial to the GSOP evaluation. The effective contribution to the GSOP is the difference between the solvation effects on neutral and oxidized molecules, eqn (11). For the TC method these values are calculated in the gas-phase. In the DC method this solvent contribution is contained in the E^{i}_{sol}(g^{i}_{sol}) energy term in eqn (1). On average the value of ΔΔG_{solv} is −1.61 eV, with a maximum value of −2.00 eV for NDI-58. This contribution does not depend largely on the method, TC or DC (ESI,† S5)
ΔΔG_{solv} = ΔG_{solv}^{+} − ΔG^{0}_{solv} | (11) |
We also evaluated the effective contribution of the thermal effects on the GSOP. Depending on the method, the thermal contribution is evaluated in solution or gas-phase for the DC or TC method respectively, for the temperature T = 298.15 K.
ΔG^{DC(TC)}_{therm}(T) = G_{therm}^{+}(g_{sol/gas}^{+},T) − G^{0}_{therm}(g^{0}_{sol/gas},T) | (12) |
Evaluation of the contribution to GSOP due to relaxation in the oxidized state, the main difference between the adiabatic and vertical approach, will follow in the next section.
ΔE^{ox} does include solvation effects (relaxation due to solvation as well) but it neglects the thermal contribution and the relaxation from the neutral to the oxidized geometry. The ΔE_{ox} data are close to the absolute experimental values with the MAD being (0.05 (RMSD = 0.18) eV). For ΔE^{ox} with R^{2} = 0.91, the strength of the linear correlation is slightly weaker than for the adiabatic approach. Comparing the differences between the two approaches, as discussed above, the thermal contributions (ESI,† S5) are rather small and consistent and should not affect the correlation strongly. In contrast, the difference coming from the relaxation due to oxidation is less consistent. On average, the energy difference due to using the oxidized or neutral geometry, both optimized in solution (gas) is 0.12 eV (0.15 eV) (ESI,† S5). The lower correlation for ΔE_{ox} is caused by a few dyes that show a larger change in energy between two oxidation states, especially NDI-5555 and NDI-555 for which the contributions due to this relaxation are 0.36 eV and 0.25 eV, respectively (ESI,† S5). However, the effect of the geometry relaxation after electron removal appears to be crucial.
Different strategies to calculate the vertical oxidation potential of the solvated dye are shown in Fig. 7. ΔE^{ox} (blue) and −ε^{DFT}_{HOMO} give very similar results and also −ε^{GW,solv}_{HOMO} gives a similar squared correlation coefficient. Compared to the other approaches, −ε^{GW,solv}_{HOMO} overestimates the GSOP considerably, while ΔE^{DFT}_{ox} (blue) and −ε^{DFT}_{HOMO} underestimate it. This depends of course on the choice of the formalism to relate the reference electrode scale to the absolute scale. However, it is known that for the gas phase, partially self-consistent GW as well as G_{0}W_{0} based on a functional with a high percentage of Hartree–Fock exchange often leads to overestimated IPs for organic molecules compared to experiment.^{99} On the other hand, using −ε^{DFT}_{HOMO} from B3LYP, due to lower percentage of Hartree–Fock exchange should rather lead to underestimation of the oxidation potential.
Fig. 7 Computed vertical GSOPs with ΔE^{ox} (blue), −ε^{DFT}_{HOMO} (red) and −ε^{GW,solv}_{HOMO}(green) compared to the experimental oxidation potential (dashed line) vs. vacuum. |
The correlation with experiment becomes considerably better when the missing thermodynamic contributions are considered. Adding the energy lowering due to geometry relaxation to −ε^{GW,solv}_{HOMO} significantly increases correlation with experiment, from R^{2} = 0.89 to R^{2} = 0.93 and for −ε^{DFT}_{HOMO} and ΔE^{ox}, the squared correlation coefficients increase from R^{2} = 0.91 to R^{2} = 0.96. These numbers should always be interpreted with the error bars in the experimental data in mind: differences in R^{2} of 0.01 (as for the difference between −ε^{DFT}_{HOMO} and −ε^{GW,solv}_{HOMO}) are negligible, but the effect of structural relaxation is significant.
In summary, the data presented in this section suggests that the vertical approximation to GSOP is not reliable. It is irrelevant what particular method is used to calculate the vertical GSOP of the solvated dye. Reliable correlation with experiment is only achieved when the relaxation of the geometry after oxidation is accounted for. This is due to the fact that this contribution can be quite large for certain systems, but almost negligible for others, as illustrated for the example of NDI-5555. However, investigating the effect of thermal contributions, we found that they are very small and that their influence is negligible. In practice, this is very important since the frequency calculations required to calculate these contributions are computationally very demanding.
As already alluded to, ΔG^{DC}_{COSMO} and ΔG^{TC}_{COSMO} are by far the computational most demanding strategies, which is due to the expensive frequency calculations. Despite their good accuracy, they are clearly not suitable for large-scale compound screening. Calculation of the G_{0}W_{0} QP energies is much cheaper but still more demanding than all other DFT-based strategies without frequency calculations. However, most of the time is spent on the QZ6P calculation. As explained in appendix 7, the GW calculations can also be performed with a TZ3P calculation only, without impairing the accuracy significantly. With the TZ3P calculation only, −ε^{GW}_{HOMO} becomes competitive with ΔE^{DFT}_{ox} and ΔG^{TC}_{COSMO-RS} and all of these strategies could be considered for large scale-screening as well. Finally, ΔG^{screening}_{COSMO-RS} is the computationally cheapest approach which is due to the fact that the geometry optimizations are performed on the GFN-xTB level. Since it provides an accuracy comparable to the involved ΔG^{DC}_{COSMO} and ΔG^{TC}_{COSMO} strategies, it is our method of choice for the screening of the dyes, described in the next subsection (Table 2).
Approach | MD | MAD | RMSD | R ^{2} |
---|---|---|---|---|
a MD stands for the mean deviation; MAD stands for the mean absolute deviation, RMSD stands for the root mean squared deviation; R^{2} is squared correlation. | ||||
ΔG^{DC}_{COSMO} | −0.28 | 0.28 | 0.30 | 0.94 |
ΔG^{TC}_{COSMO} | −0.45 | 0.45 | 0.46 | 0.95 |
ΔG^{TC}_{COSMO-RS} | −0.34 | 0.34 | 0.36 | 0.96 |
ΔE^{ox} | −0.13 | 0.15 | 0.18 | 0.91 |
−ε^{DFT}_{HOMO} | −0.05 | 0.1 | 0.13 | 0.91 |
−ε^{GW,solv}_{HOMO} | 0.20 | 0.22 | 0.24 | 0.89 |
−ε^{GW,solv,geo}_{HOMO} | 0.03 | 0.08 | 0.12 | 0.93 |
ΔG^{screening}_{COSMO-RS} | −0.34 | 0.34 | 0.37 | 0.96 |
Fig. 8 The alignment between oxidation energy levels (right axis) and potentials (left axis)of the DS-PECs components TiO_{2} CB edge (blue) and Ru-based^{100} WOC LOC (purple) with promising derivatives of four core molecules: NDI, PDI, PTI1 and PTI2 (green lines). The area under the dashed black line represents the desired GSOP values to match the highest oxidation potential of the chosen WOC. |
To better understand the effect of the root structure and the substituents, we looked at the range of GSOPs for the screened dye classes and at substituents that have the most extreme effect on the GSOP.
The unsubstituted NDI has a GSOP of 7.16 V vs. vacuum, which is lowered by substitution to a range of 6.61 V (NDI-4) to 4.40 V (NDI-5556) vs. vacuum. For the PDI, the GSOP without substitution is calculated to be 5.89 V vs. vacuum, while for the derivatives the GSOP range is from 5.64 V (PDI-4) to 4.41 V (PDI-6556) vs. vacuum. In both cases the considered substitutions are therefore lowering the GSOPs. This is not the case for the PTI1 and PTI2 core structures, where the calculated GSOPs are 6.31 V and 6.39 V respectively. For the PTI1 core the mono-substituted PTI1-4 and PTI1-1 have GSOP values of 6.49 V and 6.46 V vs. vacuum respectively, which are higher than that of the root structures. All other substituents are decreasing the GSOP value compared to the root structure, PTI1, with as lowest the PTI1-76 with a GSOP value of 5.07 V vs. vacuum. For the six PTI2 derivatives the lowest value is reached for PTI2-6, 5.49 V vs. vacuum, while the highest value is computed for monosubstituted PTI2-4, 6.03 V vs. vacuum.
For each of the root structures the ethyl alkoxide (4) substituent gives rise to the smallest GSOP shift and consequently the highest GSOP. For the NDI, PDI and PTI2 cores it decreases the value of GSOP the least compared to other substituents, while for PTI1 it increases the value of GSOP compared to unsubstituted PTI1. The substituents that appear to decrease the GSOP value the most are derived from ethanamine (5), pyrrole (6) and 4-ethynyl-N,N-dimethylaniline (7).
In Fig. 8 we show the calculated GSOP of 1340 dyes. The molecules under the dashed line fit the criterion that the GSOP is higher than the highest oxidation potential of the water oxidation catalyst, 5.65 V vs. vacuum (ESI,† S6). Out of the 118 dyes that fulfill the GSOP criterion 90 dyes are NDI derivatives, 21 PTI1 derivatives, 6 PTI2 derivatives and the final one is the PDI root structure. In contrast to this root structure, none of almost 500 PDI derivatives can fit this criterion. Of the suitable dyes the most common substituents in descending order are 4, 2 and 1. The substituents 5 and 6 appear less, while substituents 3 and 7 do not appear at all. These two substituents, derived from bithiophene and 4-ethynyl-N,N-dimethylaniline, are absent in the list of suitable dyes are decreasing the value of GSOP beyond the GSOP limit and thus prohibit the electron transfer from the WOC.
Considering only the catalyst criterion leaves a large number of candidate dyes, but one should keep in mind that the dye needs to also be able to inject an electron into the TiO_{2}, which means that a promising dye should have an ESOP lower than the TiO_{2} conduction band. The leads to an additional criterion that can be used to further narrow down the set of potentially interesting dyes.
If we approximate the ESOP as GSOP + λ_{max}, the energy corresponding to the ESOP needs to be higher than the energy corresponding to the TiO_{2} CB. As mentioned in the Introduction, for an optimal rate of the electron transfer ESOP should be approximately 0.3 V higher than the CB edge. Fig. 9 shows energy levels corresponding to GSOP and the lowest excitation energies (from ref. 17) for the dyes that fulfil absorption criteria. The purple line is the upper limit for energy that corresponds to the dyes GSOP. The dyes above the blue line satisfy the criteria that the sum GSOP + λ_{max} is higher than the energy corresponding to ESOP limit, 4.36 V vs. vacuum (TiO_{2} CB is 4.00 V vs. vacuum). A trivial way to include the expected deviations of the computed values compared to the experimental ones, is to shift the criteria for the value of MD. The dotted line, (parallel to the purple horizontal line) is the limit for GSOP including the underestimations of experimental oxidation potentials, the MD of about 0.3 V. The dashed line (parallel to the blue diagonal line) takes into account also the overestimation of the experimental absorption properties, which has a MD of about 0.4 V. The dyes in between the purple line and dashed line therefore represent the set of dyes that are suitable to this particular example system, taking into account errors in the computations. This reduces the set of suitable dyes to 66 dyes. From this set we show the dyes that have the highest oscillator strength in Table 3. The full list of these dyes is also indicated in the ESI,† S6.
Name | GSOP | λ _{max} | osc.str. |
---|---|---|---|
PTI1-5 | 6.00 | 2.90 | 0.66 |
PTI1-144 | 5.88 | 2.67 | 0.58 |
PTI1-444 | 5.81 | 2.67 | 0.56 |
PTI1-544 | 5.91 | 2.67 | 0.55 |
PTI2-4 | 6.03 | 2.70 | 0.55 |
PTI2-44 | 5.82 | 2.64 | 0.54 |
PTI1-54 | 6.02 | 2.72 | 0.53 |
NDI-125 | 5.78 | 2.57 | 0.45 |
NDI-225 | 5.84 | 2.49 | 0.44 |
PTI1-15 | 5.73 | 2.49 | 0.44 |
The dyes that have been proposed as suitable for panchromatic sensitization of the photoelectrode in DS-PECs have been further characterised for their redox properties, in particular GSOP. Using the ΔG^{screening}_{COSMO-RS} strategy, the GSOP is evaluated for 1340 dyes. For the system of choice, which is a TiO_{2} based photoelectrode and Ru-based catalyst^{100} as a WOC, there are 118 dyes that fulfil the given criteria.
Method | TZ2P | TZ3P |
---|---|---|
G _{0} W _{0}@PBE0(HF = 0.4) | 0.91 | 0.92 |
G _{0} W _{0}@LRC ω PBEH | 0.92 | 0.91 |
evGW@PBE0(HF = 0.4) | 0.87 | 0.92 |
evGW@LRCωPBEH | 0.89 | 0.90 |
qsGW_{0}@PBE | 0.88 | 0.85 |
qsGW | 0.86 | 0.91 |
For optimal agreement with experiment, the basis set error needs to be removed from the QP energies as well, using eqn (10). However, for the molecules considered here, basis set limit extrapolation increases the correlation of G_{0}W_{0}PBE0(HF = 0.4) with experiment only marginally from R^{2} = 0.92 to R^{2} = 0.93. This is due to the fact that CBS limit extrapolation results in a more or less constant shift of the QP energies, which then leaves the correlation coefficient mainly unaffected.
Footnote |
† Electronic supplementary information (ESI) available: Experimental data set (S1), R^{2} sensitivity (S2), calculated values of GSOP (S3), PDI-0000 outlier analysis (S4) and different contributions to GSOP (S5), screening results (S6). See DOI: 10.1039/d1cp04218a |
This journal is © the Owner Societies 2022 |