Why the Perfectly Symmetric Cobalt‐Pentapyridyl Loses the H2 Production Challenge: Theoretical Insight into Reaction Mechanism and Reduction Free Energies

Researchers have extensively investigated photo‐catalytic water reduction utilizing Cobalt‐based catalysts with polypyridyl ligands. While catalysts exhibiting distorted polypyridyl ligand demonstrate higher H2 production yields, those with ideal octahedral coordination display poor performance. This outcome suggests the crucial role of ligand framework in catalytic activity, yet reasons behind the disparity in H2 production rates for catalysts with octahedral geometries remain unclear. We theoretically examined the water reduction mechanism of Co‐based poly‐pyridyl catalyst, CoPy5, having perfect octahedral coordination. We clarified the effect of octahedral coordination by utilizing each intermediate step of ECEC mechanism. We determined spin states, solvent response, electronic structures, and reduction free energies. CoPy5 with perfect octahedral coordination, alongside its distorted counterparts, exhibit similar spin states as the reaction progresses through each intermediate step. However, the first reduction free energy obtained for the CoPy5 is slightly higher than that of its distorted counterparts. Following the second protonation, resulting H2 molecule experiences limited diffusion from the Co center due to the compact structure of the CoPy5, which blocks the Co center for the next H2 production cycle. Catalysts having distorted octahedral geometries facilitate fast removal of H2 into the solvent. Thus, the reaction center becomes immediately available for subsequent H2 production.


Introduction
Catalytic splitting of water into H 2 and O 2 in electrolysis requires an extrinsic energy supply.Through the agency of photocatalytic water reduction mechanism, H 2 can be also produced using solar energy. [1]Breaking down of water into H 2 and O 2 through photo-catalytic water splitting is an alternative and conspicuous method, by virtue of redefine H 2 energy as an entirely renewable energy source, thereby benefit from solar energy instead of electrical energy. [2]In order to produce molecular hydrogen, employing photo-catalytic water splitting mechanism can put into practice by means of heterogeneous and homogeneous systems.During the recent years the studies have been revolved around to comprehend the homogeneous systems and improving these systems are in high demand. [3]ccordingly, in order to achieve large scale production of H 2 along with solar radiation, usage of low costs homogeneous catalysts has draw attention in recent years. [3]s the main principle, photosensitiser at the reaction medium is excited by a photon and then transfers the electron to the water reduction catalyst.Electrons and the protons of water molecule generates H 2 molecule at the reaction center of the photo-catalyst.In the meantime the oxidized photosensitiser then receives electron from electron donor at the medium and returns to the ground state.However, it is in fact easier and applicable to implement these two processes separately. [4]The studies in literature mostly depend on exclusively water reduction or water oxidation.In order to prevent the complexity, over the course of photocatalytic H 2 production, instead of using catalyst for water oxidation, sacrificial electron donors have been used and it has been showed that donor molecules are oxidized without effecting the water reduction reaction mechanism. [5]t is been a contemporary topic to investigate H 2 production mechanism of transition metal based molecular catalysts in aqueous reaction environment. [3]Even though they increase reaction activation, it is economically not advantageous to use expensive transition metals for large scale production. [6]Instead of precious metals, e. g.Pt, Ru, in this work less expensive metal Cobalt (Co) based catalysts have been investigated.In the literature many different ligands tuneful with Co-based catalysts have been investigated, such as porphyrin, [7] cobaloxime [8] and polypyridine ligands. [4,9]Among them, by virtue of having high stability, high catalytic activity and strong ligand field, polypyridyl ligands lead the role. [10,11]It is revealed that catalytic activity can be affected by number, symmetry, and configurations of pyridyl groups around the Co center. [12]A range of results from performed experiments have shown that ligand coordinations of polypyridyl groups have profound impact on H 2 production rate. [13]ueyriaux et al. [14] studied tetra-pyridyl cobalt complex having distorted octahedral geometry where bipyridine ligands occupies the the basal plane and two water molecules occupy the axial.By achieving 443 TONs H 2 /Co their study suggests that the ECEC reaction sequence is more likely to occur with compared to the ECCE mechanism where E stands for reduction and C stands for protonation steps of the H 2 production reaction.
The cobalt based tetra-pyridyl catalyst having strongly distorted octahedral geometry has been synthesized by Alberto et al. [13] and experimental studies have been carried out to ascertain its water reduction performance.The results showed that the distorted tetra-pyridyl ligands around the Co center lead to superior H 2 turnover numbers.Bachmann et al. [15] experimentally compared the H 2 production rates of catalysts having tetra or penta pyridyl ligands which lead to distorted structures or ideal octahedral geometry around Co center, see Figure 1(d) and (e) for the structures.Experimental results showed that the catalyst having perfect octahedral geometry possesses 1200 TONs H 2 /Co in the reaction environment with 5 μM catalyst, 0.5 mM Re based photosentizer ([Re(py)(bpy)](CO) 3 ] + ), 1 M ascorbic acid/ascorbate buffer at pH 4.1 and under 385 nm LED. [15]Whereas, strongly distorted Co based catalyst having tetrapyridyl ligand exhibits 9000 TONs H 2 /Co in the reaction environment with 0.1 μM catalyst, 0.5 mM Re based photosentizer, and a 0.5/0.5 M ascorbic acid/ ascorbate buffer at the same pH and under the same LED light source. [13]As results suggest, strongly distorted Co based catalysts display significant activity compared to their perfectly symmetric counterparts.This study demonstrates that the configurations of pyridyl ligands around Co center have an effect on the catalytic activity of the system.A theoretical study carried out by Gurdal and Iannuzzi [9] compared the H 2 production reaction mechanism of Co based catalysts having penta and tetra pyridyl ligands and distorted octahedral structures.Despite experimental and theoretical studies so far it is still unclear why the ligand systems with similar electronic structures but different geometrical properties lead significantly different H 2 production rates.
The recent studies demonstrate that the information available in literature is deficient in detail mechanistic understanding of H 2 production using Co based poly-pyridyl catalysts having perfect octahedral geometry.Although experimental studies suggest higher activities linked to the distorted structures, the reason and connection between ligand structures, arrangement around the reaction center and catalytic performances have not been fully resolved.In order to fill this gap, we theoretically investigated the H 2 production mechanism of CoPy 5 catalyst, having five pyridyl ligands and perfect octahedral coordination around the Co center by employing Ab-initio Molecular Dynamics (AIMD), Density Functional Theory (DFT), as well as Free Energy Perturbation Theory.We focused on the ECEC reaction mechanism which is shown to be the probable H 2 production cycle for Co based catalyst having tetra-pyridyl ligand. [16]We modeled all intermediate steps of the reaction as well as different Co spin multiplicities for the same intermediate step in order to shed light on the possible H 2 production mechanism.The effects of arrangement of the pyridyl ligands and perfect coordination around the Co center which are responsible for having deficient H 2 production rate with respect to the distorted catalysts were discussed.

Computational Details
AIMD simulations have been performed for modeling intermediate states of the ECEC mechanisms of H 2 production through water splitting.Open source CP2K simulation package have been used in all simulations. [17]PBE density functional [18] in generalized gradient approximation (GGA) formalism was employed for the AIMD simulations.Goedecker-Teter-Hutter (GTH) potentials were applied for the estimation of core electron interactions with the valence shell and nucleus. [19]Valence electrons were modeled explicitly and valence shells of Co, N, C, O and H contain 17, 5, 4, 6 and 1 electrons, respectively.DZVP-MOLOPT basis set was used for all atomic kinds. [20]For auxiliary plane wave basis set, a cutoff of 400 Ry was utilized.Dispersion interactions were taken into consideration by applying Vydrov and Van Voorhis vdW density functional, in the revised form (rVV10). [21,22] Periodic boundary conditions and spin polarization were always applied.
For the CoPy 5 complex, AIMD simulations were carried out in a box defined as cubic with explicit water environment.The CoPy 5 catalyst was first solvated in 215 water molecules and the simulation volume was relaxed by performing AIMD simulations for approximately 20 ps in the isothermal-isobaric ensemble (NPT).Cubic simulation box volume was determined as 6163.28Å 3 .Following the determination of the simulation box size, each intermediate step were modeled by applying AIMD simulations in the canonical ensemble (NVT) for approximately 20 ps.Time step was set to 0.5 fs.Canonical sampling through velocity rescaling (CSVR) [23] thermostat with a time constant of 100 fs was applied in order to keep the simulation temperature at 300 K.All AIMD data of the reaction steps have been publicly provided in the Zenodo repository. [24]FT simulations have been applied to the selected AIMD snapshots of the intermediate steps and the electronic structure alterations of each intermediate steps were calculated using PBE0/rVV10 hybrid density functional. [25]The same approach is also used for calculating reduction free energies of the reduction steps of the ECEC mechanism by applying Free Energy Perturbation Theory.PBE0 Hartree-Fock exchange fraction parameter was set to 0.25.

Catalytic Cycle for H 2 Production
Figure 1(a) shows the heterolytic ECEC H 2 production cycle, while charges and Co spin multiplicities along the reaction pathway is depicted in Figure 1(b) where m and q represents multiplicity and total charge, respectively, as used by the m (CoPy 5 ) + q notation.The ECEC mechanism consists of subsequent electron and proton transfers to the reaction center which includes in total 2 electron and 2 proton transfers. [26] based penta-pyridyl ligand having perfect octahedral arrangement owns quartet spin multiplicity and + 2 charge at the initial state of the reaction, see Figure 1(b).Following the first electron injection, Co(I) is formed with a charge of + 1. [16,27] Accordingly, for the following step, two expectative spin states for the Co would be singlet or triplet.For the ECEC mechanism, both triplet and singlet spin states have been modeled in order to determine the most probable spin state of the catalyst that undergo the first protonation.At the ECEC mechanism, following the first protonation Co(III)-H hydride forms and at this state, the possible spin multiplicities of the Co(III) are known as quintet, triplet or singlet.Previous studies suggest that the pyridine ligands create a strong ligand field around the transition metal centers. [28]According to this, the energy levels between d orbitals of the transition metals is high at molecules which contains pyridyl ligands and in parallel with this, transition metal owns low spin multiplicity. [28,29]Consequently, in this study quintet or triplet spin states for the Co(III)-H intermediate step were neglected.Following the second electron injection, the catalyst can have doublet or quartet spin states and a charge of + 1.Both spin states have been modeled in order to determine the most appropriate spin state, depending on various traits e. g., potential energy, stability, coordination of the protons with the Co.Following the second protonation, H 2 production cycle is completed and diffusion of the H 2 to the solvent environment should be observed.

CoPy 5 Catalyst Under Consideration
CoPy 5 is composed of poly-pyridyl ligand consists of five pyridyl groups connected to each other via À OH groups and has got an perfect octahedral shape as shown in Figure 1(c).CoPy 5 provides one coordination site for water molecules to bind to the Co center.It is experimentally shown that as a water reduction catalyst, CoPy 5 , preserves its stability during the reaction. [9,15]Figure 1(d) and (e) represents the CoaPPy and CoaTPy catalysts having distorted octahedral geometry.Their H 2 production performances and water reduction reaction mechanisms have been studied experimentally and theoretically. [9,13,15]

Calculation of the Reduction Free Energies
We applied thermodynamic integration scheme and generalized Marcus Theory in order to calculate reduction free energies.More details on the applied theories and more applications can be found in Refs.[30-32].The first and the second reduction steps under consideration for calculating reduction free energies are as follows: First e À injection: Second e À injection: 1 ðCoPy 5 -HÞ þ2 þ e À ! 4 ðCoPy 5 -HÞ þ1 These expressions represent half reactions in which the number of electrons fluctuates between n in the oxidized state O and n + 1 in the reduced state R. The difference in free energy between the R and O states, denoted as ΔA RÀ O , can be calculated by the ratio of the respective partition functions.This can be further expressed as the ensemble average over the exponential of the vertical energy difference , where R N represents the instantaneous configuration of N particles.Therefore, by employing the free energy perturbation approach, [33] we can estimate the free energy difference as follows: Here, β = 1/KT, and ::: h i M denotes the ensemble average obtained from the Hamiltonian of state M.
By adopting the principles of the Marcus theory, [34,35] which suggest a linear response of the solvent to the electron transfer, we can approximate the difference in free energy by calculating the average of the ensemble averages of the vertical energy gap at both the reduced and the oxidized states as follows: Using this technique, it is possible to estimate the free energy difference from just two equilibrium simulations. [36]Under the assumption of a linear regime in the solvent's response, the probability distributions of the energy gap should exhibit a symmetric Gaussian shape and the parabolic free energy profiles possess equal curvature: [37,38] A M ðDEÞ ¼ with the curvature parameter λ being identified as the free energy of solvent reorganization: Hence, by analyzing the ensemble averages of the vertical energy gap obtained through AIMD sampling, we can determine the reduction free energy and the solvent reorganization energy.In order to evaluate the validity of the linear solvent response approximation in the electron transfer process, we also compared the variances of the energy gap distributions.
It is important to consider potential errors arising from finite size effects in charged systems when estimating solvent reorganization energy. [39]In our simulations, we have applied periodic boundary conditions and utilized the Ewald lattice summation method, [40,41] This method effectively eliminates any net charge by employing a homogeneous background charge.Previous studies have evaluated the residual errors related to solvated ions and determined them to be around 5 kJ/mol. [42]sults and Discussions We first analyze the AIMD simulations of each intermediate states of the ECEC reaction mechanism.Then, we provide electronic structure analysis of the selected snapshots of the intermediate states using hybrid density functional level of theory.Band gap fluctuations and molecular orbital distributions are analyzed for the selected AIMD snapshots.Following, calculated reduction free energies of the first and the second reduction steps are provided using free energy perturbation theory and Marcus Theory of electron transfer.

Intermediate Steps of the ECEC Reaction Mechanism
The The first step of the ECEC reaction mechanism is the first electron injection in which the system is reduced by gaining an electron.Along with the first electron injection, the total charge of the system changes as + 1 the catalyst can have either singlet or triplet spin multiplicities.We continued our simulations by injecting one electron to the system and model the first reduction step of the catalytic cycle.We simultaneously ran two AIMD-NVT simulations having triplet and singlet spin states.Figure 2(a) demonstrates a snapshot of the triplet spin state obtained after the first electron injection to the system.The simulations have been carried out around 18 ps and potential energy of the system is well preserved and temperature fluctuates at around 300 K.The O of the closest water molecule, O wat1 , preserves its distance to the Co center along the AIMD-NVT trajectory and Co-O wat1 bond dissociation have not been observed.H atoms of two different water molecules approach to the reaction center from top of the catalyst, specified as H wat2 and H wat3 .As it is depicted in Figure 2 2(c) and (d).Distance of the Co-H wat2 decreases from 4.5 Å to 3 Å.As simulation time evolves the distances again increases to 4.5 Å and decreases again periodically.H wat2 is providing a coordination with the reaction center and with the solvent by rotational moves.The atomic distance between Co-H wat2 is on the decline until 2.5 Å in conjuction with the simulation time increasing to 17 ps.The distance between Co-H wat3 , on the other hand, has large fluctuations at around 5.2 Å up to 14 ps of the simulation time.As simulation time reaches around around 18 ps, the Co-H wat3 distance becomes around 4.7 Å.As a summary, for the triplet spin state we have not observed any conformational change for the O wat1 which is very structured in the close vicinity of the reaction center.AIMD trajectory of this step is provided as a video (V.2) in the Supplementary Information (SI).Previous simulations carried out for the triplet spin state of the strongly distorted CoaTPy and CoaPPy catalysts also do not show rotation of the closely coordinated water molecule even though in the vicinity of the Co center they provide more open structures with respect to the CoPy 5 . [9]he results obtained from the singlet spin state simulations for the first electron injection step are provided at Figure 3(a).The primary significant alteration is that the closest water molecule, O wat1 , which previously provided a close coordination with the Co center changes its conformation shortly after the beginning of the simulation yielding a Co-H + bond.Along with the electron injection to the Co center, strength of the interactions between Co and O wat1 decreases and the proton of the water molecule (H wat1 ) binds to the more negative Co center.In addition to this, H of one of the other solvent water molecule, H wat4 , approaches to the reaction center.The distance between the Co-O wat1 established as 2 Å at the beginning of the simulation increases to 3.2 Å along with new Co-H wat1 bond formation, see Figure 3 sharply decreases from 4.5 Å to 2 Å due to the conformational change related to the Co-H wat1 bond formation, see Figure 3(c).AIMD trajectory of this step is provided as a video (V.1) in the SI.

(b). The distance between Co-H wat1
These results demonstrate that the singlet spin state model which is a follow up of the first electron injection leads to a bond formation between Co-H wat1 .Previous simulations carried out for the singlet spin state of the Co based distorted catalysts having tetra-pyridyl (CoaTPy) or pentapyridyl ligands (CoaPPy) showed that the water molecule which is close to Co rotates at around 5 ps and 15 ps of the simulation time, respectively. [9]For the CoPy 5 we observed such kind of rotation at around 6 ps of simulation time.As also suggested by a previous experimental and theoretical mechanistic study carried out for CoaTPy, although the triplet spin state of the Co(I) intermediate has lower potential energy than the its counterpart having singlet spin, the ECEC H 2 production reaction mechanism proceeds with the singlet spin state. [16]As a result, we proceed to the first protonation with the singlet spin state of the system.
The second intermediate step of the ECEC mechanism is the first protonation, see Figure 4(a).The proton is provided to the Co center from the closest proton of the water molecule, H wat1 , and OH À anion formed in the solvent.In order to preserve the charge of the system, an extra proton added to the simulation box away from the Co center which yields H 3 O + ion formation.A snapshot taken from the AIMD-NVT simulation of 1 (CoPy 5 ) 2 system is given in Figure 4(a).The proton binds to Co center taken from the closest water molecule and the distance between Co and the proton of the closest water molecule is depicted in Figure 4(b).In the beginning of the simulation the proton (H wat1 ) moves towards the OH À of the same water molecule which gives rise to the increased Co-H wat1 distance.However, as simulation evolves the Co center grabs the proton and Co-H wat1 distance becomes around 1.45 Å.As the simulation proceeds, the distance between oxygen atom of the OH À and its proton (H wat1 ) increases from 1 Å to around 2.5 Å which also implies that the bond between OH À and the proton of the water molecule dissociates, see Figure 4(c).At the same time at around 2.5 ps, the proton of the closest water molecule establishes a bond with the Co center of the water reduction catalyst, see Figure 4(b).Besides, the remained OH À forms hydrogen bonding interactions with the surrounding water molecules and diffuses to the solvent, we refer V.3 in the SI for the corresponding AIMD trajectory.The proton, on the other hand, binds to the Co center and remains stable until the end of the simulation.As the results suggest, the proton is supplied to the Co center by the closest water molecule and Co(III)-H hydride is stable.
The next step of the ECEC mechanism is the second electron injection and this intermediate step can proceed either with doublet or quartet spin states.Within our study, we modeled both doublet and quartet spin states for the second electron injection step, see Figure 1(b).At Figure 5(a) distances between the Co-H + and H wat1 as well as Co-H + and H wat2 are plotted along the simulation time for the doublet spin state.Note that here H + represents the proton binded to the Co and H wat1 and H wat2 are the water hydrogens closest to the proton binded Co center.In the beginning of the simulation, the distance between Co-H + and H wat1 was around 4 Å.As the simulation continues, the mentioned distance decreases to 1.7 Å.The distance between Co-H + and H wat2 , on the other hand, increases from 2.4 Å to 3.5 Å. Although, the simulation is continued to 17 ps, the alteration of H wat2 and H wat1 coordination with the Co(III)-H hydride seems to continue.The AIMD-NVT snapshot provided from the trajectory of the second electron injection with the doublet spin state is given at Figure 5(b).From the given results, it can be inferred that for doublet multiplicity state only one proton from the solvent can coordinate with the Co(III)-H hydride.While H wat1 is coordinating with the Co-H + , the other H wat2 is interacting with the solvent.As the same vein, while H wat2 is coordinating with the Co-H + , H wat1 is interacting with the solvent through hydrogen bonding.
The results of the quartet spin state simulations are depicted in Figure 5(c) which shows the distances between Co-H + and H of the closest water molecule, Co-H wat1 , as well as  Co-H + and the the second closest water molecule, Co-H wat2 , Figure 5(d) shows the corresponding.The Co-H wat1 distance first increases from 1.5 Å to 2.5 Å and then decreases from 2.5 Å to around 1.5 Å again and remains stable until the end of the simulation.Simultaneously, the Co-H wat2 distance decreases from 3.5 Å to 1.5 Å.These results give a conclusion that both protons of water molecules can provide a coordination with the Co center, as opposed to the doublet spin state.Due to increased possibility of available protons around the Co center, the H 2 production reaction mechanism is likely to continue with quartet spin multiplicity for the second protonation.
The fourth and the last step of the ECEC mechanism is the second protonation.This step gives rise to production of the H 2 molecule and the multiplicity as well as the charge of the system should, in principle, revert back to the initial state as, 4 (CoPy 5 ) + 2 -H 2 .Figure 6(a) shows the bond distances between protons of the H 2 molecule (in purple color) which is formed with the second protonation step.Ideally, with the second protonation step the two protons should bind to each other and move away from the Co center by diffusing towards the solvent.Conformably to the ideal situation, we observed that the distances between two protons fluctuates at around 0.8 Å and preserves its stability until the end of the simulation, which signifies that the H 2 molecule is, indeed, formed.However, due to the presence of the surrounding water molecules and geometric nature of the catalyst, H 2 molecule struggles with moving away from the Co center and the distance between Co center and the H 2 molecule, indicated in orange color in Figure 6(a), remains stable at around 1.9 Å.At Figure 6(b) a snapshot is provided from the 9 ps of the simulation time.It is, indeed, interesting that we could not observe diffusion of the H 2 from the vicinity of the Co center, see V.4 in the SI.Previous simulations carried out for the strongly distorted CoaTPy and CoaPPy catalysts, on the other hand, already show diffusion of the H 2 molecules through the solvent shortly after the beginning of the simulations. [9]ur results suggest that the perfect octahedral coordination of the CoPy 5 catalyst and its more closed structure around the Co center with respect to the distorted CoaTpy and CoaPPy have an adverse effect on the diffusion of the produced H 2 molecules.The H 2 production through water splitting using CoPy 5 could be hindered by the diffusion of the products from the reaction center which explains less activity observed for the perfectly octahedral CoPy 5 with respect to the distorted structures.

Structural Changes Observed for the Catalyst
Along the simulations of the H 2 production reaction, catalyst undergoes structural changes and accordingly these changes effect the distances between the pyridyl groups and the Co center.Alterations of mentioned distances has an important place to specify the geometry of catalyst.As provided in the Table S1 in the SI, distances between pyridine rings and Co center slightly decreases following the first reduction step, which is more effective for the singlet spin state.Depending on reducing the system, the catalyst becomes slightly more compact.The first protonation has no significant effect on the bond lengths between Ns and Co.Following the second electron injection, a significant elongation of bonds between Ns and Co center were observed for both doublet and quartet spin multiplicities, however the effect is more visible for the quartet spin state.Following the second protonation step of the ECEC mechanism, the production cycle comes to an end, charge and multiplicity of the catalyst return to the initial state, and catalyst restores its initial structure.

Electronic Structure Analysis of Each Intermediate Steps
Atomic rearrangements and solvent response influence the electronic structures, thus in order to understand the electronic structure responses, we observed the evolution of electronic structures for each intermediate state in the reaction cycle.We calculated band gaps and the highest occupied (HOMO) and the lowest unoccupied molecular orbital (LUMO) distributions of each intermediate step by applying PBE0-rVV10 level of  theory.For each intermediate step 60 snapshots were taken from the corresponding AIMD trajectories and HOMO and LUMO distributions were visualized using VMD. [43]HOMO and LUMO evolution of the snapshots were provided as videos in the SI.
In Figure 7 HOMO and LUMO distributions of the selected snapshots of each intermediate steps were provided.For the initial state, 4 (CoPy 5 ) + 2 , HOMO and LUMO representations were obtained from 16 ps of the AIMD-NVT trajectory, see Figure 7(b).HOMO is mainly located on the Co center, Ns, and the closest water molecule, while LUMO is distributed mostly on pyridyl rings belonging to N APy , N RPy1 and N RPy2 .Following the first electron injection, the distribution of the HOMO starts from N APy including Co center, proceeds at axial direction and ends at where the water molecule binds to Co. LUMO, on the other hand, is distributed mainly on Co center and pyridyls which are located at the axial position, see Figure 7(c).For the first protonation, we observe large HOMO-LUMO energy difference.
For instance energy gaps at 15.4 and 18.8 ps are calculated as 3.68 and 1 eV, respectively.HOMO corresponding to large band gap state at 15.4 ps is located on right hand side pyridyl rings as well as on hydroxy group, see Figure 7(d).While for the small band gap state, obtained at 18.8 ps, HOMO is distributed to the solvent.Corresponding LUMO states are distributed to the solvent, as shown as videos in the SI.As shown in the Figure 7(e), at 15 ps HOMO is distributed on the solvent for the second electron injection step which has a relatively small band gap of 1.02 eV.However, at 18 ps, HOMO is mainly located on the catalyst, specifically on the Co, N APy , and N RPy2 atoms as well as on the the closest water molecule.

Reduction Free Energy
H 2 production through the ECEC reaction mechanism has two reduction steps.As soon as the first electron is injected to the system, the system can have ether singlet or triplet spin states.As solvent rearrangement suggests, reaction proceeds to the first protonation step with the singlet spin multiplicity.Therefore, for the first electron injection step, reduction free energy calculations were carried out only for the singlet spin state.While, for the first electron reduction, the oxidized limit is the initial state of the system, 4 (CoPy 5 ) + 2 , following the reduction, the reduced limit is the 1 (CoPy 5 ) + 1 .Following the second electron injection, 4 (CoPy 5 ) + 1 is the reduced state where 1 (CoPy 5 ) + 2 is the oxidized state.
We selected more than 180 snapshots from each of the AIMD-NVT trajectories and calculated the vertical energy differences by applying PBE0 together with rVV10 vdW density functional in order to increase accuracy of the calculated reduction free energies.We analyzed the energy differences between oxidized and reduced states and corresponding energy differences, ΔE, over the production time in picoseconds for the first reduction step as well as the second reduction steps of the ECEC mechanisms, see Figure 8 As depicted in the Figure 8(a), the ensemble averages of the vertical energy differences DE h i R and DE h i O for the singlet spin state are estimated as À 4.79 eV and 0.79 eV, respectively.Using Eq. ( 2), the reduction free energy, ΔA RÀ O , for the singlet spin state of the first electron injection was calculated as, À 2 eV.In the case of distorted structures, on the other hand, reduction free energies of the first reduction reaction have been calculated as À 1.81 eV and À 1.95 eV for the singlet spin state of the CoaPPy and CoaTPy, respectively. [9]Thus, the reduction free energy associated with the CoPy 5 is slightly higher than that of its distorted counterparts.
The reorganization free energies were also calculated for each electron injection step in order to determine the energy needed for solvent molecules to reorganize themselves when any electronic effects occur at the system.Reorganization energy of singlet spin state of the first electron injection is estimated as À 2.79 eV by using Eq. ( 4).This large value can be explained by the water molecules situated around the catalyst

Figure 1 .
Figure 1.(a) Visual representations of the ECEC reaction mechanism for the CoPy 5 .Co(II) depicts initial state of the catalyst and by reduction and protonation intermediate steps H 2 is formed.(b) Intermediate steps modeled in this work.m and q respresents multiplicity and total charge, respectively, as used by the m (CoPy 5 ) + q notation.(c) Ball and stick representation of the CoPy 5 catalyst modeled in this work.N atoms in the pyridyl rings are labeled as N APy , N RPy1 , N RPy2 , N LPy1 , and N LPy2 .(d) CoaPPy and (e) CoaTPy catalysts having distorted octahedral arrangements and studied in Refs.[9,13,15].Color code: green: C, white: H, purple: Co, blue: N, red: O.
initial state of the CoPy 5 owns quartet spin multiplicity and the system has + 2 total charge.AIMD simulations in the canonical ensemble (NVT) have been carried out to simulate initial state of the catalyst in explicit water environment.Water molecules disperse homogeneously around the water reduction catalyst.Oxygen of one of the water molecule binds to the Co center which is labeled as O wat1 , see Figure 2(a).The distance between the Co center and the O wat1 fluctuates at around 2.1 Å throughout the AIMD-NVT simulations of the initial state.The distances between N atoms within the pyridyl ligands of the CoPy 5 and Co center fluctuate at around 2.1 Å and during the AIMD simulations bond dissociation between the ligands and the Co center have not been observed.

Figure 2 .
Figure 2. (a) Snapshot taken from the AIMD-NVT simulations of the triplet spin state 3 (CoPy 5 ) + 1 at 17 ps of simulation time.(b) Alterations along the simulation for distances between O of the closest three water molecules and Co center of the catalyst.(c) Fluctuations of distances between Co center and H of the second closest water molecule (H wat2 ) as well as.(d) Co center and H of the third closest water molecule (H wat3 ).

Figure 3 .
Figure 3. (a) Snapshot taken from the AIMD-NVT simulations of the singlet spin state 1 (CoPy 5 ) + 1 at 19.63 ps of simulation time.(b) Alterations along the simulation for distances between O of the closest two water molecules and Co center of the catalyst.(c) Fluctuations between Co center and H atoms of the two closest water molecules.

Figure 4 .
Figure 4. (a) A snapshot taken from the AIMD-NVT simulation of the first protonation step, 1 (CoPy 5 ) + 2 , at 15 ps of simulation time is given.(b) After the first protonation, 1 (CoPy 5 ) + 2 , the atomic distance between the Co and the bonded proton provided by the nearest water molecule (H wat1 ) is plotted in blue, Co-H wat1 .(c) The atomic distance between the OH À anion and the proton bonded to the Co is plotted in green, O wat1 -H wat1 .

Figure 5 .
Figure 5. (a) Co-H wat1 and Co-H wat2 distances along the AIMD for 2 (CoPy 5 ) + 1 which are depicted by purple and orange, respectively, where Co represents the Co-H + bond.(b) A snapshot illustrating the 2 (CoPy 5 ) + 1 state at 16 ps of simulation time, (c) Co-H wat1 and Co-H wat2 distances along the AIMD for the 4 (CoPy 5 ) + 1 are depicted by blue and red, respectively, where Co represents Co-H + bond.(d) A snapshot depicts the 4 (CoPy 5 ) + 1 state taken at 17 ps of simulation time.

Figure 6 .
Figure 6.(a) H + -H + distances and Co-H 2 distances along the AIMD trajectory are plotted in purple and orange colors, respectively.(b) A snapshot is illustrated following the second protonation step.
(a) and (c), respectively.The frequencies of energy differences are depicted in Figure 8(b) and (d).

Figure 8 .
Figure 8.(a) The vertical energy differences versus production time (ps) shown for the first electron injection.(b) Frequencies of the vertical energy differences of the first reduction reaction, red shows oxidation and blue shows reduction reactions.(c) The vertical energy differences versus production time (ps) of the second electron injection, (d) Frequencies versus the vertical energy differences for the second electron injection step.In (a) and (c), turquoise lines represents the energy differences of the reduced limit and accumulative average shown with dark blue, whereas pink lines represents energy differences of oxidized limits and its accumulative averages shown with red.
(b), the closest water molecule, O wat1 , preserves its coordinate during the simulations and its distance to the Co center which fluctuates at around 2.1 Å.There beside, distances between Co and the second closest water molecule, Co-O wat2 , and the third closest water molecule, Co-O wat3 fluctuate widely at around 4 Å.The coordinations between H atoms of the corresponding water molecules, H wat2 and H wat3 , and the Co center of the catalyst are depicted at Figure