Two‐Channel Model for Electron Transfer in a Dye‐Catalyst‐Dye Supramolecular Complex for Photocatalytic Water Splitting

Abstract To improve the performance of dye‐sensitized photoelectrochemical cell (DS‐PEC) devices for splitting water, the tailoring of the photocatalytic four‐photon water oxidation half‐reaction represents a principle challenge of fundamental significance. In this study, a Ru‐based water oxidation catalyst (WOC) covalently bound to two 2,6‐diethoxy‐1,4,5,8‐diimide‐naphthalene (NDI) dye functionalities provides comparable driving forces and channels for electron transfer. Constrained ab initio molecular dynamics simulations are performed to investigate the photocatalytic cycle of this two‐channel model for photocatalytic water splitting. The introduction of a second light‐harvesting dye in the Ru‐based dye‐WOC‐dye supramolecular complex enables two separate parallel electron‐transfer channels, leading to a five‐step catalytic cycle with three intermediates and two doubly oxidized states. The total spin S=1 is conserved during the catalytic process and the system with opposite spin on the oxidized NDI proceeds from the Ru=O intermediate to the final Ru‐O2 intermediate with a triplet molecular 3O2 ligand that is eventually released into the environment. The in‐depth insight into the proposed photocatalytic cycle of the two‐channel model provides a strategy for the development of novel high‐efficiency supramolecular complexes for DS‐PEC devices with buildup and conservation of spin multiplicity along the reaction coordinate as a design principle.


S1. Expanded reaction pathways
Scheme S1. Expanded Reaction Pathways of the Two-Channel Model by a Ru-based dye-WOC-dye System a a Schematic Structure of the starting intermediate 1 ([(cy)Ru II bpy(H2O)] 2+ −(NDI)2) complex (indicated shortly as 1 (NDI1-[Ru II -OH2] 2+ −NDI2)) on the top of the scheme, as explicitly shown in the inset. It is assumed that each light flash induces an electron injection (golden arrows) from the NDI1/NDI2 to the semiconductor electrode or to the next stage in a tandem cell, leading to the photooxidation of NDI1/NDI2: NDI1/NDI2 → NDI1 +• /NDI2 +• . All the possible pathways for the dye−WOC−dye complex after introducing an extra NDI dye are considered in this scheme. In the current work we only focus on the case involving the cophotooxidation of two NDI dyes and the stable intermediates involved as shown in Scheme 1 are indicated in red. The transient states indicated in black are not investigated in this work.
To obtain a realistic description of the catalytic reaction step, the solvent was explicitly the AIMD simulations to get accurate predictions of the catalytic reaction, which was carried out with the CPMD program. [5] Two reaction coordinates were considered in our constrained MD simulations, corresponding to the distances between Hiii/Hiv and the oxygen atom of its neighboring water molecule (see the atomic labeling in Scheme S1 and Figure S1). The reaction coordinates d(Hiii←O) as well as d(Hiv←O) are constrained to a series of fixed values x in range of 1.4 − 1.05 Å simultaneously after the initial equilibrium simulation and subsequent photooxidation of two NDI dyes. In this way, we assume that one attacking water molecule approaches to Hiii and another attacking water molecule to Hiv at the same time. The photooxidation of the NDI dyes (NDI1 and NDI2) was mimicked by removing two electrons from the simulation box after the initial equilibration simulation of the dye−WOC−dye system, after which a free MD (FMD) simulation of around 1.1 ps at room temperature was performed to equilibrate the oxidized state 3 (NDI1 +• −[Ru II −OH2] 2+ −NDI2 +• ) with the total spin S = 1. According to the results of our simulations, one unpaired α electron (↑) is observed to localize on NDI1 and one unpaired α electron (↑) on NDI2 in the system (see inset (i) in Figure S1). This spin density confirms that Spin density on WOC

S5
The shortening of d(Hiii←O) and d(Hiv←O) from 1.4 Å to 1.1 Å induces the electron transfer from the WOC to the oxidized NDI dyes (see Figure S1a) The Oi−Hiii and Oi−Hiv bonds (see inset (iii) in Figure S1 for the atomic labeling) finally break when we further shorten the Hiii•••Oiv and Hiv•••Ov distances to 1.05 Å (see Figure S1b and S1c), which occurs almost at the same time as the accomplishment of the electron transfer (see Figure S1a). No back-transfer of either an electron or a proton is observed after the release of the constraints from the system following the constrained 1.05 Å simulation. This confirms the stability of the final product 3 (NDI1−[Ru IV =O] 2+ −NDI2) (S = 1) with two unpaired α electrons localized on the WOC after the first and second catalytic PCET steps (see inset (ii) in Figure S1). All these results indicate that the first and second catalytic steps are able to proceed and complete at the same stage after the photooxidation of two NDI dyes.  Figure S2. Table S1 summarizes the key thermodynamic parameters extracted from the free energy profiles  The obtained free energy profiles for these two PCET processes reported in Figure S5 show similar activation free energy barriers ∆G* ≈ 3.9 kcal mol -1 (0.17 eV) and ∆G* ≈ 4.6 kcal mol -1 (0.20 eV), corresponding to the reaction coordinates d(Hiii←O) and d(Hiv←O) respectively in the two-channel model (see Table S1). This is consistent with the comparable activation barriers of the first and second catalytic PCET steps in the one-channel model [2,7] . The maximum of the free energy profile corresponds to a reaction coordinate d(Hiii←O)/d(Hiv←O) of 1.14/1.14 Å,

S3.1 Geometry Optimization at DFT Level
The OPBE exchange-correlation functional [1] and the TZP (triple-ζ polarized) Slater-type basis set [2] were employed in the geometry optimization of the initial state of the dye−WOC−dye complex. The OPBE functional has shown to be accurate in describing transition-metal complexes, including Ru-based WOCs. [8] In the geometry optimization, the continuous solvation model (COSMO [4] ) for water was used. These calculations are performed with the ADF software package. [3]

S3.2 Simulation Box
To obtain a realistic description of the catalytic reaction step, the solvent was explicitly introduced in the simulations. The solvent environment for the CPMD simulations was generated using Discovery Studio 2.5. [9] The solvent was equilibrated for 0.2 ns using the TIP3P model implemented in the CHARMM force field and CFF partial charge parameters at 300 K, [10] while the dye−[WOC] 2+ −dye complex was kept fixed. The volume was then adjusted using constant pressure for 0.2 ns, after which the system was further allowed to evolve with constant volume for 2 ns. Periodic boundary conditions are applied with a time step of δt = 5 a.u. (1 a.u. = 0.0242 fs).

S3.3 Free Energy Profile
To estimate the activation free energy barrier of the catalytic reaction step involving the O-O bond formation that is unlikely to occur spontaneously during the typical AIMD simulation time scale, constrained MD and the so-called Blue Moon approach were employed as a rare event simulation technique. [6] The reaction coordinate (in this case the distance between two oxygen atoms Oi and Oii, d(Oi←Oii), as shown in Scheme 1) is constrained to a series of fixed values x in range of 2.5 − 1.6 Å after the initial equilibrium simulation and subsequent photooxidation of two NDI dyes along this reaction pathway. A time-averaged constraint force <λ>x for each value of the reaction coordinate x is obtained, which should be equal to zero at an equilibrium or transition state. Based also on our previous work on a similar supramolecular complex, we can safely assume that for a value of the reaction coordinate d(Oi←Oii) = 3.0 Å the system is in an equilibrium state with <λ>3.0 Å equals to zero. [7] The activation free energy barrier for this catalytic step is then established by interpolating the mean forces with a 100-S9 point Akima splines function and integrating the signed forces <λ>x along the reaction path. [11] Trajectory analysis and visualization for the CPMD output were carried out using VMD. [12] S3.4 Reaction Rate The computed activation free energy barrier can be used to evaluate to what extent the geometry modification accelerates the rate of the third water oxidation step involving the O−O bond formation. According to transition state theory [13] , the reaction rate (k) determined by the activation energy barrier (∆G*) can be expressed as Where ∆G* represents the activation free energy barrier, kB, h, R and T are the Boltzmann constant, the Planck constant, the universal gas constant and thermodynamic temperature, respectively. One should keep in mind that in the DFT-based MD simulations protons are treated classically and thus proton tunneling effects are neglected. In the current calculation, only the activation energy barrier is considered as a main factor governing the reaction rate. S10 S4. Molecular orbital and electronic structure of 3 (NDI1 +• -[Ru IV =O] 2+ -NDI2 +• )  Table S3).  Figure S5. Spin density integrated over the part of the simulation box including the WOC along the free MD trajectory after the photooxidation of two NDI dyes. Inset shows the spin density isosurface computed at a snapshot taken at ~0.12 ps, clearly indicating two unpaired α electrons (↑ in green) localized on the catalyst, one unpaired β electron on NDI1 (↓ in purple), and one unpaired α electron on NDI2 (↑ in green  Figure S6. Time evolution of the geometrical parameters (blue line) and (black line) along the constrained MD trajectory with d(Oi←Oii) = 1.8 Å. The inset shows the schematic structure of the first two water molecules along the hydrogen-bonding network coordinated to the oxygen ligand. The time range is consistent with that in Figure 1. S14 S8. Time-averaged constraint forces Figure S7. Time-averaged constraint force represented by the Lagrangian multiplier <λ> computed for each constrained MD simulation as a function of the reaction coordinate d(Oi←Oii) in the two-channel model. The Akima splines (100 points) is used to interpolate the mean forces. The final intermediates corresponding to the MD simulations for the one-channel and two-channel models are both indicated. The time-averaged constraint forces obtained in the one-channel model from a previous study is also presented for comparison (see Ref. [7]). The error bars indicate the standard deviations.