Open Access Article
This Open Access Article is licensed under a
Creative Commons Attribution 3.0 Unported Licence

Automated assessment of redox potentials for dyes in dye-sensitized photoelectrochemical cells

Jelena Belić a, Arno Förster a, Jan Paul Menzel b, Francesco Buda b and Lucas Visscher *a
aDepartment of Chemistry and Pharmaceutical Sciences, Vrije Universiteit Amsterdam, De Boelelaan 1083, 1081 HV, Amsterdam, The Netherlands. E-mail:
bLeiden Institute of Chemistry, Leiden University, Einsteinweg 55, P.O. Box 9502, 2300 RA, Leiden, The Netherlands

Received 15th September 2021 , Accepted 30th November 2021

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.

1 Introduction

Hydrogen is seen as a symbol of sustainable energy source, even though its production currently mostly relies on fossil fuels.1,2 To make production of this energy carrier more sustainable,3–5 a promising route is to split water with a dye-sensitized photoelectrochemical cell (DS-PEC) driven by solar energy.6,7 Such developments capitalize on the advantages of dye-sensitization, in particular the possibility to tune the dye's properties by small structural adjustments.7

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 TiO2). 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.

image file: d1cp04218a-f1.tif
Fig. 1 The alignment between redox potentials of the DS-PEC's components: conduction band (CB) edge of the anode (blue), ground state oxidation potential (GSOP) and excited state oxidation potential (ESOP) of the dye (green) and highest oxidation potential (HOP) of the water oxidation catalyst (purple).

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 devices7,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 anode10–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 injection18 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 generalized22,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 energies33 and is therefore increasingly used for energy level alignment in photo-catalytic interfaces34,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+ + eanode) 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 E1/2 of the reaction is, at standard conditions, a very good estimate for the standard electrode potential: E1/2E°,40 with E° being the driving force for electrochemical work.

As we are primarily interested in the GSOP here, the measured E1/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.

2 Methods

In this section we will discuss the possible computational strategies for GSOP determination and the experimental data that can be used to validate them. We will start with the two adiabatic strategies directed towards full determination of the Gibbs free energy, consider in the second section approaches that use the vertical ionization potential, and end with the selection of the experimental reference data.

2.1 GSOP – adiabatic approach

In the adiabatic approach, the GSOP is defined as the absolute difference between the solution-phase Gibbs free energies of the products and reactants, Dye and Dye+ + e. It describes the energy change between the molecular species in their respective solvation and thermal equilibrium conditions. The Gibbs free energy of the electron depends on the choice of the statistical mechanical formalism – which does not matter as long as it is consistent with the formalism of the reference electrode, which we will address later. For the Fermi–Dirac statistics this value is −3.62 kJ mol−1 (−0.0375 eV).42,43 Solution-phase Gibbs free energies Gisol(gisol), where the index i labels the neutral (0) or oxidized (+) forms of the dye and g stands for the optimized geometry of the species, can be calculated using DFT. We thereby distinguish between the dye's energy Eisol(gisol) and thermal contributions Githerm(gisol,T) at the temperature T of interest:
Gisol(gisol) = Eisol(gisol) + Githerm(gisol,T)(1)
The subscripts hereby indicate the phase of the species: solvated (sol) or in gas phase.

The thermal contribution is the most demanding term in the calculation as it requires determination of the vibrational frequencies of the molecule. The Githerm(gisol,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 model44 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 (ΔGoxsol) 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 COSMO45 as CS model, ΔGDCCOSMO is calculated using eqn (2):

ΔGDCCOSMO = Gsol+(gsol+) + Ggas(e) − G0sol(g0sol)(2)
Gsol+(gsol+) and G0sol(g0sol) are solution-phase Gibbs free energies of the oxidised and neutral molecular species, respectively.

image file: d1cp04218a-f2.tif
Fig. 2 The thermodynamic cycle. The molecule's total Gibbs free energies are denoted with G and molecule's geometry with g in brackets; the subscript is the phase and in superscript is the oxidation state of the species. Solvation processes are the vertical and oxidation processes are the horizontal reactions. Next to the arrows are changes in Gibbs free energies for each process. The free electron is by electrochemical convention always in gas phase.

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 Gigas(gigas) (following eqn (1) for gas phase), the Gibbs free energy of solvation computed at the gas phase geometry gigas, and the energy difference between the solvated and gas phase structure:

Gisol(gigas) = Gigas(gigas) + ΔGisolv(gigas) + [Egas(gisol) − Egas(gigas)](3)

The equation for the adiabatic GSOP computed with the TC method and COSMO, ΔGTCCOSMO is then similar to eqn (2) and reads:

ΔGTCCOSMO = Gsol+(ggas+) + Ggas(e) − G0sol(g0gas)(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 robust48 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. Ho47 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,

Gisol,CRS(gigas) = Eigas(gigas) + ΔGiCRS,solv(gigas)(5)
Inserting the solution-phase Gibbs free energy defined in eqn (5) into eqn (4), the adiabatic GSOP using the TC pathway in combination with COSMO-RS is obtained.

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 ΔGscreeningCOSMO-RS.

2.2 GSOP – vertical approach

The Gibbs free energy difference should correspond exactly to the GSOP measured in the CV experiments as the time scale of these experiments is long enough to allow for full relaxation. Computing the GSOP in the vertical approximation and neglecting vibrational effects is computationally attractive and can be considered a valid approximation if the system is not in thermodynamical equilibrium, which will be the case when the electron transfer process is fast enough. Physically, the relaxation of the neutral geometry after electron removal is then neglected and only relaxation of the electrons as a response to the removal of an electron to the system is accounted for ref. 32. In this regime one can consider three different strategies.
2.2.1 Calculation as total electronic energy difference. The first approach is to compute total electronic energy for both the neutral and oxidized species at the optimized geometry of neutral molecule in solvent. The total energy difference
ΔEox = Esol+(g0sol) − E0sol(g0sol)(6)
is then taken as approximation to the GSOP.
2.2.2 Calculation as KS-DFT HOMO energy. The simplest approach would be to identify the negative of the HOMO energy, −εHOMO, from a KS-DFT calculation with the vertical ionisation energy in exact KS-DFT. While for any finite system −εHOMO exactly equals the ionisation energy in exact KS-DFT,58 errors obtained from common approximations to the exact exchange–correlation potential vxc can be of the order of several eV.29,30,59 This failure is a direct consequence of the incorrect long-range behaviour of the vxc(r) (exponentially instead of 1/|r|) in the local density approximation (LDA) and generalized gradient approximations (GGA).28 On the other hand, the exact exchange (eex) potential shows the correct long-range behaviour60 and hybrid functionals, combining GGAs with a fraction of eex give typically much better ionisation energies. Despite these obvious shortcomings, the KS-DFT approach is computationally very efficient since it only requires to perform a single geometry optimization in solution.
2.2.3 Calculation as GW HOMO energy. In a single-particle energy framework, one can also go beyond DFT and calculate a vertical approximation to GSOP using MBPT. Central to MBPT is Dyson's equation61
image file: d1cp04218a-t1.tif(7)
where εi is a KS eigenvalue, ϕi denotes a molecular orbital, ωi is the exact one-electron energy and Σ is the so-called self-energy. In practice, the GWA to the self-energy31 is often used. In the GWA, Σ is calculated as the convolution (in frequency space) of the interacting Green's function G and the dynamically screened interaction W, which is related to the unscreened Coulomb interaction vc by
W0(ω) = [εRPA(ω)]−1vc, εRPA = 1 + vcπ0(ω),(8)
where π0 is the independent particle polarizability in the RPA. Eqn (7) is usually simplified by evaluating Σ with a non-interacting Green's function, G0, instead of the interacting one. In the G0W0 approximation62 Σ is additionally approximated as diagonal, so that the eqn (7) reduces to a set of independent non-linear equations,
Σii(ωi) = ωiεi,(9)
with Σii = ϕi|Σ(ω)|ϕi. Thus, the QP energies are obtained as a perturbative correction to the KS eigenvalues. In eigenvalue-only self-consistent GW (evGW), the QP energies are additionally updated until self-consistency is reached.62 Yet another option is to map Σ to a static, non-local and Hermitian exchange–correlation potential,63,64 which then defines a non-linear eigenproblem, much like in KS-DFT, with the only difference that the potential is a functional of G0. This approach is referred to as QP self-consistent GW (qsGW). The GW approximation can be implemented with almost quadratic scaling using localized atomic orbitals,65–67 which makes at least G0W0 competitive with hybrid functional calculations in terms of computational cost in the same basis set.67 Self-consistent variants are more expensive, but the calculations typically converge in 5–10 iterations.68 However, GW calculations converge much slower to the complete basis set (CBS) limit.33,69,70 For accurate QP energies, it is usually necessary to perform calculations using correlation consistent basis sets of at least triple- and quadruple-ζ quality and to extrapolate to the CBS limit.69–71

2.3 Experimental data set

To validate the methods discussed above we compared the calculated GSOP values with experimental values for E1/2. The set used for validation is a subset of all considered dye cores (Fig. 3) for which we only selected those GSOP measurements that were done under similar experimental conditions. The substituted naphthalene-diimides (NDI),72–75 perylene-diimide (PDI)76 core and substituted PDI cores76,77 with different number and different type of substituents are used for validation of theoretical predictions.
image file: d1cp04218a-f3.tif
Fig. 3 Structures of the unsubstituted cores considered for screening.

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 Fc/Fc+ 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).

image file: d1cp04218a-f4.tif
Fig. 4 Structures of the substituents attached to the cores.

To quantify the predicting ability of the methods we will use mean deviations and squared correlation – R2. 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 R2 – 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 R2 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 R2 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 R2 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.

3 Computational details

All computations were performed with the Amsterdam Modeling Suite (AMS) from Software for Chemistry and Materials (SCM),78 versions 2019.302 and 2020.1, except for the GW calculations which were performed with a modified development version of 2020.1. Initial geometries for the core and the library of substituents were prepared with the SCM Graphical User Interface (GUI). The long alkyl chains were replaced with ethyl groups. Generation of NDI and PDI derivative structures was done with the Compound Attachment Tool (CAT).79 CAT is a Python code employing, among others, the Python Library for Automating Molecular Simulation (PLAMS)80 for generating structures. The initial molecular structures were pre-optimized with the DFTB381 method using the 3-ob parameter set.82–85 Convergence to a stable minimum of the nuclear potential energy surface was checked by calculating only the lowest frequency normal modes.86 We allowed for some numerical noise and considered structures with imaginary frequencies above −20 cm−1 for further characterisation. The workflow for the strategies above in general consists of geometry optimization (GO) in some cases followed by frequencies (normal modes calculation) and single point (SP) calculation(s). SP calculations were done with DFT or GW. Depending on the strategy, solvation effects were included in the geometry optimization calculation or in the following single point calculation (ESI, S3). The DFT calculations (GO and SP) were done using the B3LYP functional,87 with the TZ2P basis set.88

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 COSMO90 or COSMO-RS91 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.30292 settings were used: the cavity construction is Delley type,93 atomic radii are the corresponding van der Waals radii from the MM3 method by Allinger94 (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 (G0W0) and partially self-consistent GW: G0W0@PBE0(HF = 0.4) (G0W0 performed on top of a PBE095,96 calculations with 40% exact exchange), G0W0@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 (qsGW0@PBE). All GW calculations were performed using the TZ2P and TZ3P71 basis sets. For G0W0@PBE(HF = 0.4) we additionally used the QZ6P basis set71 and extrapolated the IP to the CBS limit using

image file: d1cp04218a-t2.tif(10)
where εQZn (εTZn) denotes the value of the QP energy using QZ6P (TZ3P) and NQZbas and NTZbas denote the respective numbers of basis functions (in spherical harmonics so that there are 5d and 7f functions).71

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 qsGW0 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 G0W0 and evGW we set NumericalQuality Good and for the qsGW0 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 G0W0 and evGW calculations and Dependency Bas = 5 × 10−3 for qsGW0.

4 Results

4.1 Adiabatic GSOP

We compare the strategies for obtaining the Gibbs free energy of the oxidation reaction with two pathways: TC and DC with COSMO. For the TC strategy we additionally consider the COSMO-RS method. The performance of these methods is determined for 13 out of 14 molecules from the initial data set. We excluded the PDI-0000 molecule as it turned out to heavily influence the correlation with the experimental reference data and thus presents a clear outlier. The unexpectedly large errors seen for this molecule may be caused by uncertainties in selecting its most appropriate conformation98 or by its sensitivity to CS model parameters (ESI, S4). Both aspects make this molecule unsuited for getting accurate results with the automated procedure for initial structure generation and subsequent calculations that we currently employ. Table 1 shows that calculation timings are similar for the ΔGTCCOSMO and ΔGDCCOSMO method. Both methods perform equally well for the given set of experimental results, Fig. 5. The TC method has a slightly higher squared correlation coefficient (R2 = 0.95) than the DC methods (R2 = 0.94). The values from the TC method have an almost consistent shift from the corresponding experimental values, with a mean absolute deviation of 0.45 eV (RMSD = 0.46 eV). The DC method is a bit closer to the absolute experimental values with a MAD of 0.28 eV (RMSD = 0.30 eV). The ΔGTCCOSMO-RS absolute deviation of 0.34 eV (RMSD = 0.36 eV) is as large as for other two methods, while the squared correlation is higher than both methods employing only COSMO, R2 = 0.96. The ΔGscreeningCOSMO-RS, where the geometry optimisation is performed using the SQM method (GFN1-xTB), performs as well as ΔGTCCOSMO-RS. The correlation between these two methods is almost 1 and for both methods the correlation with experiment is 0.96.
Table 1 The strategies considered in this publication together with required calculations, different methods, and corresponding computation times. Method performance was accessed for DFTc with B3LYP functional in combination with other methods. Timings denote elapsed times measured in CPU hours and are for unsubstituted PDI with 40 atoms
Approach Calculationsa 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).
ΔGDCCOSMO GO(n,s) + f DFT 390 1066
GO(o,s) + f 676
ΔGTCCOSMO 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
ΔGTCCOSMO-RS GO(n,v) DFT 9.6 32.7
+SP(n,s) 1.4
GO(o,v) 19.2
+SP(o,s) 2.5
ΔEDFTox GO(n,s) DFT 9.6 12.1
+SP(o,s) 2.5
ε DFTHOMO 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
ΔGscreeningCOSMO-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

image file: d1cp04218a-f5.tif
Fig. 5 The correlation of adiabatic GSOP computed with ΔGTCCOSMO (pink), ΔGDCCOSMO (grey) and ΔGTCCOSMO-RS (purple) methods to the experimental oxidation potential vs. vacuum (dashed line).

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 ΔGTCCOSMO 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 Eisol(gisol) energy term in eqn (1). On average the value of ΔΔGsolv 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)

ΔΔGsolv = ΔGsolv+ − ΔG0solv(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.

ΔGDC(TC)therm(T) = Gtherm+(gsol/gas+,T) − G0therm(g0sol/gas,T)(12)
The effective ΔGDC(TC)therm(T) contribution for DC and TC, has an average of 0.04 eV and 0.03 eV respectively (for the absolute value difference). So, frequencies calculated in the same environment are very close in value for neutral and oxidized dyes and their contribution cancel out. Considering the cost of frequency calculations listed in Table 1, for GO(n,s) + f and just GO(n,s), the frequencies are 97% of GO(n,s) + f time. We therefore conclude that by neglecting this small thermal contribution, workflows for this type of dyes can be sped up significantly.

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.

4.2 Vertical approximations to GSOP

The alternative to the full adiabatic GSOP is to approximate it by a vertical oxidation potential. Apart from experimental considerations, such an approximation can always be justified if the effect of the geometry relaxation as well as other thermal contributions are negligible. The validity of this approximation can be assessed directly by comparing the Gibbs free energy for the oxidation reaction obtained with DC, ΔGDCCOSMO to the electronic energy difference ΔEox: the calculations are performed in the same manner and we can focus on the physical approximations. The comparison is shown in Fig. 6.
image file: d1cp04218a-f6.tif
Fig. 6 The GSOP computed as Gibbs free energy of the oxidation using the DC pathway (pink) and the vertical approach employing total electronic energy difference (blue) compared to the experimental oxidation potential (dashed line) vs. vacuum.

ΔEox 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 ΔEox data are close to the absolute experimental values with the MAD being (0.05 (RMSD = 0.18) eV). For ΔEox with R2 = 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 ΔEox 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. ΔEox (blue) and −εDFTHOMO give very similar results and also −εGW,solvHOMO gives a similar squared correlation coefficient. Compared to the other approaches, −εGW,solvHOMO overestimates the GSOP considerably, while ΔEDFTox (blue) and −εDFTHOMO 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 G0W0 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 −εDFTHOMO from B3LYP, due to lower percentage of Hartree–Fock exchange should rather lead to underestimation of the oxidation potential.

image file: d1cp04218a-f7.tif
Fig. 7 Computed vertical GSOPs with ΔEox (blue), −εDFTHOMO (red) and −εGW,solvHOMO(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,solvHOMO significantly increases correlation with experiment, from R2 = 0.89 to R2 = 0.93 and for −εDFTHOMO and ΔEox, the squared correlation coefficients increase from R2 = 0.91 to R2 = 0.96. These numbers should always be interpreted with the error bars in the experimental data in mind: differences in R2 of 0.01 (as for the difference between −εDFTHOMO and −εGW,solvHOMO) 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.

4.3 Timings

Since we want to find a method which is suitable for large-scale screening for potentially thousands of compounds, the CPU times of the different approaches are an important consideration. We use the unsubstituted PDI with 40 atoms as a representative example. Some of the systems considered here are about twice as large, but our conclusions do not change for these systems. All calculations presented in this subsection have been performed on a single 2.2 GHz intel Xeon (E5-2650 v4) node (broadwell architecture) with 24 cores and 128 GB RAM. Timings are shown in Table 1.

As already alluded to, ΔGDCCOSMO and ΔGTCCOSMO 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 G0W0 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, −εGWHOMO becomes competitive with ΔEDFTox and ΔGTCCOSMO-RS and all of these strategies could be considered for large scale-screening as well. Finally, ΔGscreeningCOSMO-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 ΔGDCCOSMO and ΔGTCCOSMO strategies, it is our method of choice for the screening of the dyes, described in the next subsection (Table 2).

Table 2 Statistical analysisa of the considered strategies compared with cyclic voltammetry measurements in dichloromethane
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; R2 is squared correlation.
ΔGDCCOSMO −0.28 0.28 0.30 0.94
ΔGTCCOSMO −0.45 0.45 0.46 0.95
ΔGTCCOSMO-RS −0.34 0.34 0.36 0.96
ΔEox −0.13 0.15 0.18 0.91
εDFTHOMO −0.05 0.1 0.13 0.91
εGW,solvHOMO 0.20 0.22 0.24 0.89
εGW,solv,geoHOMO 0.03 0.08 0.12 0.93
ΔGscreeningCOSMO-RS −0.34 0.34 0.37 0.96

4.4 Screening

We have employed the ΔGscreeningCOSMO-RS method to calculate the GSOP potentials of large number of dyes. We have applied it to set of dyes that we determined previously17 to have desirable optical properties. Among almost 2500 NDIs, PDIs, PTI1s and PTI2s derivatives around 1400 dyes fulfilled the criteria to be used as a photosensitizers in DS-PECs. All these dyes have an intense, lowest transition in the desired energy range, between the 1.35 eV which is the minimal thermodynamic requirement for water oxidation and 3.20 eV which is roughly the boundary between the visible and the UV energy range. During the water oxidation cycle the WOC goes through different oxidation states. For the GSOP of the dye the most important criteria is that it should lie higher than the WOC's highest oxidation potential (Fig. 1) to assure a favourable potential gradient for electron transfer. We thereby assume 0.1 V potential difference to be sufficient. As a trial WOC we will take the Ru-based catalyst designed by Duan et al.100 for which the most demanding catalytic step occurs at the potential 1.27 V vs. NHE (5.55 V vs. vacuum). Therefore, the lower limit for the dye's oxidation potential is at 5.65 V vs. vacuum (dashed black line in Fig. 8).
image file: d1cp04218a-f8.tif
Fig. 8 The alignment between oxidation energy levels (right axis) and potentials (left axis)of the DS-PECs components TiO2 CB edge (blue) and Ru-based100 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 TiO2, which means that a promising dye should have an ESOP lower than the TiO2 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 TiO2 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 (TiO2 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.

image file: d1cp04218a-f9.tif
Fig. 9 The lowest excitation energy and energy levels corresponding to GSOP potentials for the 1340 dyes. The color of the dot corresponds to its absorption energy and its opacity corresponds to oscillator strength normalised to 0 to 1 scale. The horizontal purple line represents the experimental GSOP limit and dotted line (parallel to purple line) includes computational deviations. The blue line represents the experimental ESOP limit (the sum of GSOP and λmax needs to exceed the TiO2 CB). The dashed line (parallel to blue line) includes computational deviation. The optimal set of dyes are in between the purple and dashed line.
Table 3 GOSP [V vs. vacuum], lowest most intense excitation λmax [eV] and oscillation strength (osc.str) values for 10 suitable dyes with the highest oscillation strengths17
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

5 Conclusions

We have computed the solution-phase Gibbs free energy, with adiabatic and vertical approaches, employing different methods. First we have built the experimental data set from the CV measurements for the derivatives of the NDI and PDI dyes to validate the computational results. Calculations of solution-phase Gibbs free energy are done at the DFT level using different CSM pathways, the DC and TC pathways, and using different CSM models – COSMO and COSMO-RS. We find that, to calculate the ground state oxidation potential for these dyes, both pathways using the COSMO model perform well. The TC path shows a slightly higher value of squared correlation with the experiment, but a higher MAD value as well. The TC pathway with the COSMO-RS model has the highest correlation with the experiment with a MAD lower than the TC pathway with COSMO. Comparison to different vertical approximation, ΔEox and one electron energy strategies, the −εDFTHOMO and −εGWHOMO, shows the importance of taking into account geometry relaxation after oxidation, as well as the inclusion of electronic solvent effects. However, other thermal effects do not play an important role when it comes to this set of molecules. Therefore, the adiabatic approach appears more suitable for screening purposes, but one should be aware that this procedure might be unfavourable for the molecules with large number of conformational isomers that exist in a small energy range and that are highly affected by solvent, such as PDI-0000. Replacing DFT geometry optimization by SQM optimization has not significantly affected the correlation with experiment for the adiabatic approach combined with the COSMO-RS model, while it reduced the cost of the strategy by a factor of eight. Therefore, the ΔGscreeningCOSMO-RS is found to be a suitable strategy for screening on a desired GSOP range for derivatives of the NDI and PDI dyes.

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 ΔGscreeningCOSMO-RS strategy, the GSOP is evaluated for 1340 dyes. For the system of choice, which is a TiO2 based photoelectrode and Ru-based catalyst100 as a WOC, there are 118 dyes that fulfil the given criteria.

Conflicts of interest

There are no conflicts to declare.


A Effect of solvent and structural relaxation on GW

The output of the GW calculations are QP energies. They correspond to the vertical electron addition and removal energies in vacuum and can not be compared to the GSOP directly. This is shown here for G0W0@PBE0(HF = 0.4) in Fig. 10, the same is observed for all other tested GW approaches. Identifying −εGWHOMO systematically overestimates the GSOP by several eV on average. Also the squared correlation is with only R2= 0.87 worse than for ΔG. Accounting for solvent effects removes most of the overestimation and improves R2 to 0.89. Taking into account the structural relaxation after electron removal improves the correlation to R2 = 0.93, which is now competitive with ΔG. Note, that this does not include other thermal contributions. As discussed above, they are very small and can be neglected. Our results are consistent with earlier work by Umari et al.37 who compared G0W0@BLYP vertical oxidation potentials to experimental data from differential pulse voltammetry in acetonitrile.101 They observed that their GW approach overestimated the experimental values by nearly one eV on average, even though it is known, that G0W0 based on a GGA starting point systematically underestimates the gas phase oxidation potential of organic molecules.99
image file: d1cp04218a-f10.tif
Fig. 10 The GSOP calculated with the GW approaches successively including more physical effects (solvation effects and geometry relaxation due to oxidation): −εGWHOMO, −εGW,solvHOMO and −εGW,solv,geoHOMO (lightest to darkest shade of green) compared to the experimental oxidation potential (dashed line) vs. vacuum.

B Comparison and validation of GW methods

Taking into account these considerations, we can compare the different GW methods to experiment. We tested all approaches using the TZ2P and TZ3P basis sets and calculated the R2 values with respect to experiment. The result of our comparison is shown in Table 4. The correlation of G0W0 was found to be more or less independent of the starting point. Eigenvalues-only self-consistency does not improve the results compared to G0W0 for TZ3P and worsens them for TZ2P. The same is observed for qsGW. qsGW0 gives much worse correlation coefficient on the TZ3P level (R2 = 0.85) than the other methods. For this reason, using one of the more expensive partially self-consistent approaches is not justified and we used the G0W0PBE0(HF = 0.4) method for comparison to the other computational methods in the main body of the text. Our conclusion that self-consistency is not needed here is of course not valid in general and in other cases it might be possible that updates in the wave-function are needed. Reordering of orbitals during the GW calculation can be an indicator for this and in such cases it might be necessary to use some kind of self-consistency in GW.36
Table 4 The squared correlation coefficients R2 compared to experiment for different GW approaches for different basis sets. Solvent corrections and the effects of geometry relaxation are included
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
qsGW0@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 G0W0PBE0(HF = 0.4) with experiment only marginally from R2 = 0.92 to R2 = 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.


This research has been financially supported by the NWO Solar to Products program (project number 733.000.007). AF acknowledges funding from the Netherlands Organisation for Scientific Research (NWO) in the framework of the Innovation Fund for Chemistry and from the Ministry of Economic Affairs in the framework of the TKI/PPS-Toeslagregeling. We acknowledge the use of supercomputer facilities at SURFsara sponsored by NWO Physical Sciences, with financial support from The Netherlands Organization for Scientific Research (NWO).


  1. The Future of Hydrogen: Seizing today's opportunities, IEA Technical Report June, 2019.
  2. Hydrogen Roadmap Europe: A sustainable pathway for the European Energy Transition, Fuel cells and hydrogen joint undertaking technical report, 2019.
  3. K. Fujishima and A. Honda, Nature, 1972, 238, 37–38 CrossRef.
  4. C. Acar and I. Dincer, Chem. Eng. Sci., 2019, 197, 74–86 CrossRef CAS.
  5. C. Decavoli, C. L. Boldrini, N. Manfredi and A. Abbotto, Rend. Lincei, 2019, 30, 469–483 CrossRef.
  6. W. Song, Z. Chen, M. Kyle Brennaman, J. J. Concepcion, A. O. T. Patrocinio, N. Y. Murakami Iha and T. J. Meyer, Pure Appl. Chem., 2011, 83, 749–768 CAS.
  7. C. Decavoli, C. L. Boldrini, N. Manfredi and A. Abbotto, Eur. J. Inorg. Chem., 2020, 978–999 CrossRef CAS.
  8. R. Katoh, A. Furube, T. Yoshihara, K. Hara, G. Fujihashi, S. Takano, S. Murata, H. Arakawa and M. Tachiya, J. Phys. Chem. B, 2004, 108, 4818–4822 CrossRef CAS.
  9. C. Li and H. Wonneberger, Adv. Mater., 2012, 24, 613–636 CrossRef CAS.
  10. J. T. Kirner and R. G. Finke, J. Mater. Chem. A, 2017, 5, 19560–19592 RSC.
  11. J. T. Kirner, J. J. Stracke, B. A. Gregg and R. G. Finke, ACS Appl. Mater. Interfaces, 2014, 6, 13367–13377 CrossRef CAS.
  12. F. Ronconi, Z. Syrgiannis, A. Bonasera, M. Prato, R. Argazzi, S. Caramori, V. Cristino and C. A. Bignozzi, J. Am. Chem. Soc., 2015, 137, 4630–4633 CrossRef CAS.
  13. M. Bonchio, Z. Syrgiannis, M. Burian, N. Marino, E. Pizzolato, K. Dirian, F. Rigodanza, G. A. Volpato, G. La Ganga, N. Demitri, S. Berardi, H. Amenitsch, D. M. Guldi, S. Caramori, C. A. Bignozzi, A. Sartorel and M. Prato, Nat. Chem., 2019, 11, 146–153 CrossRef CAS.
  14. C. E. Creissen, J. Warnan, D. Antón-García, Y. Farré, F. Odobel and E. Reisner, ACS Catal., 2019, 9, 9530–9538 CrossRef CAS.
  15. B. Van Den Bosch, J. A. Rombouts, R. V. Orru, J. N. Reek and R. J. Detz, ChemCatChem, 2016, 8, 1392–1398 CrossRef CAS.
  16. H. C. Chen, R. M. Williams, J. N. Reek and A. M. Brouwer, Chem. – Eur. J., 2016, 22, 5489–5493 CrossRef CAS PubMed.
  17. J. Belić, B. van Beek, J. P. Menzel, F. Buda and L. Visscher, J. Phys. Chem. A, 2020, 124, 6380–6388 CrossRef PubMed.
  18. J. P. Menzel, A. Papadopoulos, J. Belić, H. J. de Groot, L. Visscher and F. Buda, J. Phys. Chem. C, 2020, 124, 27965–27976 CrossRef CAS.
  19. Y. Shao, H. J. de Groot and F. Buda, ChemSusChem, 2021, 14, 479–486 CrossRef CAS PubMed.
  20. Y. Shao, J. M. De Ruiter, H. J. De Groot and F. Buda, J. Phys. Chem. C, 2019, 123, 21403–21414 CrossRef CAS.
  21. F. De Angelis, S. Fantacci and A. Selloni, Nanotechnology, 2008, 19, 424002 CrossRef PubMed.
  22. A. Seidl, A. Görling, P. Vogl, J. Majewski and M. Levy, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 53, 3764–3774 CrossRef CAS PubMed.
  23. A. Görling and M. Levy, J. Chem. Phys., 1997, 106, 2675–2680 CrossRef.
  24. W. Kohn and L. J. Sham, Phys. Rev., 1965, 140, A1133 CrossRef.
  25. P. Hohenberg and W. Kohn, Phys. Rev., 1964, 136, 864–871 CrossRef.
  26. R. E. Skyner, J. L. McDonagh, C. R. Groom, T. Van Mourik and J. B. Mitchell, Phys. Chem. Chem. Phys., 2015, 17, 6174–6191 RSC.
  27. R. W. Godby, M. Schlüter and L. J. Sham, Phys. Rev. B: Condens. Matter Mater. Phys., 1988, 37, 10159–10175 CrossRef.
  28. E. Engel, J. A. Chevary, L. D. Macdonald and S. H. Vosko, Zeitschrift für Phys. D Atoms, Mol. Clust., 1992, 23, 7–14 CrossRef CAS.
  29. E. J. Baerends, O. V. Gritsenko and R. Van Meer, Phys. Chem. Chem. Phys., 2013, 15, 16408–16425 RSC.
  30. E. J. Baerends, Phys. Chem. Chem. Phys., 2017, 19, 15639–15656 RSC.
  31. L. Hedin, Phys. Rev., 1965, 139, A796–A823 CrossRef.
  32. G. Onida, I. Nazionale, R. T. Vergata, R. Scientifica and I. Roma, Rev. Mod. Phys., 2002, 74, 601 CrossRef CAS.
  33. D. Golze, M. Dvorak and P. Rinke, Front. Chem., 2019, 7, 1–66 CrossRef.
  34. A. Migani, D. J. Mowbray, A. Iacomino, J. Zhao, H. Petek and A. Rubio, J. Am. Chem. Soc., 2013, 135, 11429–11432 CrossRef CAS.
  35. A. Migani, D. J. Mowbray, J. Zhao, H. Petek and A. Rubio, J. Chem. Theory Comput., 2014, 10, 2103–2113 CrossRef CAS.
  36. N. Marom, J. E. Moussa, X. Ren, A. Tkatchenko and J. R. Chelikowsky, Phys. Rev. B: Condens. Matter Mater. Phys., 2011, 84, 1–15 Search PubMed.
  37. P. Umari, L. Giacomazzi, F. De Angelis, M. Pastore and S. Baroni, J. Chem. Phys., 2013, 139, 0147091–0147098 CrossRef.
  38. N. Marom, T. Körzdörfer, X. Ren, A. Tkatchenko and J. R. Chelikowsky, J. Phys. Chem. Lett., 2014, 5, 2395–2401 CrossRef CAS.
  39. D. J. Mowbray and A. Migani, J. Phys. Chem. C, 2015, 119, 19634–19641 CrossRef CAS.
  40. A. J. Bard and L. R. Faulkner, ELECTROCHEMICAL METHODS Fundamentals and Applications, John Wiley & Sons, Inc, New York, 2nd edn, 2001 Search PubMed.
  41. J. Preat, C. Michaux, D. Jacquemin and E. A. Perpète, J. Phys. Chem. C, 2009, 113, 16821–16833 CrossRef CAS.
  42. J. E. Bartmess, J. Phys. Chem., 1994, 98, 6420–6424 CrossRef CAS.
  43. A. V. Marenich, J. Ho, M. L. Coote, C. J. Cramer and D. G. Truhlar, Phys. Chem. Chem. Phys., 2014, 16, 15068–15106 RSC.
  44. S. Grimme, Chem. – Eur. J., 2012, 18, 9955–9964 CrossRef CAS.
  45. A. Klamt and G. Schüürmann, J. Chem. Soc., Perkin Trans. 2, 1993, 799–805 RSC.
  46. R. F. Ribeiro, A. V. Marenich, C. J. Cramer and D. G. Truhlar, J. Phys. Chem. B, 2011, 115, 14556–14562 CrossRef CAS PubMed.
  47. J. Ho, Phys. Chem. Chem. Phys., 2015, 17, 2859–2868 RSC.
  48. J. Ho, M. L. Coote, C. J. Cramer and D. G. Truhlar, in Org. Electrochem., ed. O. Hammerich and B. Speiser, CRC Press, Boca Raton, FL, 5th edn, 2016, ch. 6, pp. 229–259 Search PubMed.
  49. C. J. Cramer and D. G. Truhlar, Acc. Chem. Res., 2008, 41, 760–768 CrossRef CAS PubMed.
  50. A. Klamt, B. Mennucci, J. Tomasi, V. Barone, C. Curutchet, M. Orozco and F. J. Luque, Acc. Chem. Res., 2009, 42, 489–492 CrossRef CAS PubMed.
  51. J. Ho and M. L. Coote, Theor. Chem. Acc., 2009, 125, 3–21 Search PubMed.
  52. J. Ho, A. Klamt and M. L. Coote, J. Phys. Chem. A, 2010, 114, 13442–13444 CrossRef CAS PubMed.
  53. A. Klamt, J. Phys. Chem., 1995, 99, 2224–2235 CrossRef CAS.
  54. A. Klamt, V. Jonas, T. Bürger and J. C. Lohrenz, J. Phys. Chem. A, 1998, 102, 5074–5085 CrossRef CAS.
  55. A. Klamt, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2011, 1, 699–709 CAS.
  56. H. Neugebauer, F. Bohle, M. Bursch, A. Hansen and S. Grimme, J. Phys. Chem. A, 2020, 124, 7166–7176 CrossRef CAS.
  57. J. Menzel, M. Kloppenburg, J. Belić, H. de Groot, L. Visscher and F. Buda, J. Comput. Chem., 2021, 42, 1–10 CrossRef.
  58. C. O. Almbladh and U. Von Barth, Phys. Rev. B: Condens. Matter Mater. Phys., 1985, 31, 3231–3244 CrossRef CAS.
  59. E. J. Baerends, O. V. Gritsenko and R. Van Meer, Phys. Chem. Chem. Phys., 2013, 15, 16408–16425 RSC.
  60. J. D. Talman and W. F. Shadwick, Phys. Rev. A: At., Mol., Opt. Phys., 1976, 14, 36–40 CrossRef CAS.
  61. F. J. Dyson, Phys. Rev., 1949, 75, 1736–1755 CrossRef.
  62. M. S. Hybertsen and S. G. Louie, Phys. Rev. B: Condens. Matter Mater. Phys., 1986, 34, 5390–5413 CrossRef CAS.
  63. M. Van Schilfgaarde, T. Kotani and S. Faleev, Phys. Rev. Lett., 2006, 96, 1–4 Search PubMed.
  64. T. Kotani, M. Van Schilfgaarde and S. V. Faleev, Phys. Rev. B: Condens. Matter Mater. Phys., 2007, 76, 1–24 CrossRef.
  65. J. Wilhelm, D. Golze, L. Talirz, J. Hutter and C. A. Pignedoli, J. Phys. Chem. Lett., 2018, 9, 306–312 CrossRef CAS.
  66. J. Wilhelm, P. Seewald and D. Golze, J. Chem. Theory Comput., 2021, 17, 1662–1677 CrossRef CAS PubMed.
  67. A. Förster and L. Visscher, J. Chem. Theory Comput., 2020, 16, 7381–7399 CrossRef PubMed.
  68. A. Förster and L. Visscher, Front. Chem., 2021, 9, 736591 CrossRef.
  69. M. J. Van Setten, F. Caruso, S. Sharifzadeh, X. Ren, M. Scheffler, F. Liu, J. Lischner, L. Lin, J. R. Deslippe, S. G. Louie, C. Yang, F. Weigend, J. B. Neaton, F. Evers and P. Rinke, J. Chem. Theory Comput., 2015, 11, 5665–5687 CrossRef CAS.
  70. A. Stuke, C. Kunkel, D. Golze, M. Todorović, J. T. Margraf, K. Reuter, P. Rinke and H. Oberhofer, Sci. Data, 2020, 7, 1–11 CrossRef.
  71. A. Förster and L. Visscher, J. Chem. Theory Comput., 2021, 17, 5080–5097 CrossRef.
  72. C. Röger and F. Würthner, J. Org. Chem., 2007, 72, 8070–8075 CrossRef.
  73. B. A. Jones, A. Facchetti, M. R. Wasielewski and T. J. Marks, J. Am. Chem. Soc., 2007, 129, 15259–15278 CrossRef CAS PubMed.
  74. S. Bhosale, A. L. Sisson, N. Sakai and S. Matile, Org. Biomol. Chem., 2006, 4, 3031–3039 RSC.
  75. R. S. K. Kishore, O. Kel, N. Banerji, D. Emery, G. Bollot, J. Mareda, A. Gomez-Casado, P. Jonkheijm, J. Huskens, P. Maroni, M. Borkovec, E. Vauthey, N. Sakai and S. Matile, J. Am. Chem. Soc., 2009, 131, 11106–11116 CrossRef CAS PubMed.
  76. P. Leowanawat, A. Nowak-Król and F. Würthner, Org. Chem. Front., 2016, 3, 537–544 RSC.
  77. M. J. Ahrens, M. J. Tauber and M. R. Wasielewski, J. Org. Chem., 2006, 71, 2107–2114 CrossRef CAS.
  78. SCM,
  79. B. van Beek and J. Belić, CAT version 0.8.7, DOI:10.5281/zenodo.3832787, 2020.
  80. PLAMS,;
  81. M. Gaus, Q. Cui and M. Elstner, J. Chem. Theory Comput., 2011, 7, 931–948 CrossRef CAS PubMed.
  82. M. Kubillus, T. Kubař, M. Gaus, J. Řezáč and M. Elstner, J. Chem. Theory Comput., 2015, 11, 332–342 CrossRef CAS PubMed.
  83. X. Lu, M. Gaus, M. Elstner and Q. Cui, J. Phys. Chem. B, 2015, 119, 1062–1082 CrossRef CAS PubMed.
  84. M. Gaus, A. Goez and M. Elstner, J. Chem. Theory Comput., 2013, 9, 338–354 CrossRef CAS.
  85. M. Gaus, X. Lu, M. Elstner and Q. Cui, J. Chem. Theory Comput., 2014, 10, 1518–1537 CrossRef CAS.
  86. P. Deglmann and F. Furche, J. Chem. Phys., 2002, 117, 9535–9538 CrossRef CAS.
  87. P. J. Stephens, F. J. Devlin, C. F. Chabalowski and M. J. Frisch, J. Phys. Chem., 1994, 98, 11623–11627 CrossRef CAS.
  88. E. Van Lenthe and E. J. Baerends, J. Comput. Chem., 2003, 24, 1142–1156 CrossRef CAS.
  89. S. Grimme, C. Bannwarth and P. Shushkov, J. Chem. Theory Comput., 2017, 13, 1989–2009 CrossRef CAS PubMed.
  90. C. Pye and T. Ziegler, Theor. Chem. Acc., 1999, 101, 396 Search PubMed.
  91. C. Pye, T. Ziegler, E. van Lenthe and J. Louwen, Can. J. Chem., 2009, 87, 790 CrossRef CAS.
  92. ADF 2019.3,
  93. B. Delley, Mol. Simul., 2006, 32, 117–123 CrossRef CAS.
  94. N. L. Allinger, X. Zhou and J. Bergsma, J. Mol. Struct.: THEOCHEM, 1994, 312, 69–83 CrossRef.
  95. C. Adamo and V. Barone, J. Chem. Phys., 1999, 110, 6158–6170 CrossRef CAS.
  96. M. Ernzerhof and G. E. Scuseria, J. Chem. Phys., 1999, 110, 5029–5036 CrossRef CAS.
  97. M. A. Rohrdanz, K. M. Martins and J. M. Herbert, J. Chem. Phys., 2009, 130, 054112-1-8 CrossRef.
  98. F. Würthner, Chem. Commun., 2004, 1564–1579 RSC.
  99. J. W. Knight, X. Wang, L. Gallandi, O. Dolgounitcheva, X. Ren, J. V. Ortiz, P. Rinke, T. Körzdörfer and N. Marom, J. Chem. Theory Comput., 2016, 12, 615–626 CrossRef CAS PubMed.
  100. L. Duan, F. Bozoglian, S. Mandal, B. Stewart, T. Privalov, A. Llobet and L. Sun, Nat. Chem., 2012, 4, 418–423 CrossRef CAS PubMed.
  101. D. P. Hagberg, T. Marinado, K. M. Karlsson, K. Nonomura, P. Qin, G. Boschloo, T. Brinck, A. Hagfeldt and L. Sun, J. Org. Chem., 2007, 72, 9550–9556 CrossRef CAS.


Electronic supplementary information (ESI) available: Experimental data set (S1), R2 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