Tuning the Proton‐Coupled Electron‐Transfer Rate by Ligand Modification in Catalyst–Dye Supramolecular Complexes for Photocatalytic Water Splitting

Abstract In view of the considerably high activation energy barrier of the O−O bond formation photocatalytic step in water oxidation, it is essential to understand if and how nonadiabatic factors can accelerate the proton‐coupled electron transfer (PCET) rate in this process to find rational design strategies facilitating this step. Herein, constrained ab initio molecular dynamics simulations are performed to investigate this rate‐limiting step in a series of catalyst‐dye supramolecular complexes functionalized with different alkyl groups on the catalyst component. These structural modifications lead to tunable thermodynamic driving forces, PCET rates, and vibronic coupling with specific resonant torsional modes. These results reveal that such resonant coupling between electronic and nuclear motions contributes to crossing catalytic barriers in PCET reactions by enabling semiclassical coherent conversion of a reactant into a product. Our results provide insight on how to engineer efficient catalyst‐dye supramolecular complexes by functionalization with steric substituents for high‐performance dye‐sensitized photoelectrochemical cells.


S1.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 and final states of WOC−dye complexes L0 -L3. The OPBE functional has shown to be accurate in describing transition-metal complexes, including Ru-based WOCs. [3] In the geometry optimization, the continuous solvation model (COSMO [4] ) for water was used. These calculations are performed with the Amsterdam Density Functional (ADF) software package. [5]

S1.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 Car−Parrinello MD (CPMD) simulations was generated using Discovery Studio 2.5. [6] 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, [7] while the [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).

S1.3 Effect of periodic boundary conditions
Periodic boundary conditions (pbc) are applied in our simulations. In plane wave based AIMD simulations the periodic boundary conditions introduce a spurious Coulomb interaction for charged systems due to the image charges. The effect of pbc for charged systems can be important when dealing with isolated molecules/clusters in the simulation box. However, because of the quite large simulation box (25.1 × 17.7 × 14.4 Å 3 ) used and the fact that the MD simulation box contains 162 water molecules that will strongly screen the spurious Coulomb interaction. The spurious effect of periodic charges is estimated to be rather small (∼0.01 eV). We can therefore conclude that the error introduced by the pbc does not affect significantly the conclusions of our work.

S1.4 Free energy profile
To estimate the activation free energy barrier of the catalytic reaction step involving the O-O bond formation that are unlikely to occur spontaneously during the typical ab initio molecular dynamics (AIMD) simulation time scale, constrained MD and the so-called Blue-Moon approach were employed as a rare event simulation technique. [8] 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 NDI 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. Considering what we learned from our previous work that the reaction coordinate d(Oi←Oii) = 3.0 Å and 1.325 Å are quite close to the initial/final equilibrium state and far from the transition state, we therefore assume that the modification of ligand R has minor effect on the position close to the initial/final equilibrium state and all the complexes L0 -L3 share the same value for the constraint forces at the reaction coordinate d(Oi←Oii) = 3.0 Å, 2.7 Å, and 1.325 Å. [9] The activation free energy barrier for this catalytic step is then established by interpolating the mean forces with a 100-point Akima splines function and integrating the signed forces <λ>x along the reaction path. [10] Trajectory analysis and visualization for the CPMD output were carried out using VMD. [11]

S1.5 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 the transition state theory [12] , the reaction rate (k) determined by the activation energy barrier (∆G*) can be expressed as (1) 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. Table S1. Calculated geometrical parameters for complexes L0 -L3.    along the free MD trajectories after the photooxidation of NDI dye. The insets show the spin density isosurface computed at a snapshot taken at the end of each free MD simulation of complexes L0 -L3, respectively, in the doublet state with two unpaired α electrons localized on the catalyst (green spin density isosurface) and one unpaired β electron on the NDI dye (purple spin density isosurface). An integrated spin density value of 1 corresponds to one unpaired β electron (↓).   The distance between Ru and H3O + , defined as an oxygen atom with 3 H within a radius of 1.2 Å, illustrating the proton diffusion during the constrained 1.8 Å MD simulations for complexes L0 -L3. The analysis of the trajectories shows that only one oxygen is in the H3O + form at any time, and the excess proton associates primarily to 2 -4 different oxygens (indicated with different colours) during the simulation. This figure clearly shows that the proton diffuses from the first solvation shell of the ruthenium center to the second or even third and fourth solvation shell rapidly within this relatively short MD timescale of ~0.6 ps. The time range is consistent with that in Figure 3. Figure S7. Snapshot taken at the end of the constrained 1.8 Å simulation for complex L2. The attacking water molecule and the neighboring water molecules forming the hydrogen-bonded chain are represented with ball & stick. The dashed blue lines indicate the hydrogen bonds. This figure clearly shows that the proton Hi has been totally released by the attacking water molecule and diffuses rapidly into the solvent bulk via a "chain" of hydrogen-bonded water molecules even within the MD timescale of ~0.6 ps.   Figure S12. Spin density integrated over half of the simulation box including the NDI dye of complex L3 along the constrained 2.5 Å MD trajectory with fixed dihedral angle θ = 91° after the photooxidation of NDI dye (blue line). Before this constrained MD simulation, a simulation of ∼0.36 ps with fixed dihedral angle θ = 91° was performed to equilibrate the solvated system. The data without constraint on the dihedral angle θ of L3 (grey line) is also presented for comparison, which is extracted from Figure 3.   a The geometry of all the complexes L0 -L3 together with the attacking water molecule are firstly optimized with the ADF program using OPBE functional and the TZP basis set. The TD-DFT calculations are then performed at the same lavel. The continuum solvation model (COSMO) is used to describe the water environment. The distance between the two oxygen atoms Oi and Oii (d(Oi-Oii), Å) is fixed to certain values which is taken around the transition state according to Figure 4. The first excitation is mainly related to the transition from HOMO (SOMO WOC) to LUMO (SOMO dye). SOMO represents the singly occupied molecular orbital.