UvA-DARE (Digital Academic Unraveling the mechanism of biomimetic hydrogen fuel production - a first principles molecular dynamics study

The Fe 2 (bdt)(CO) 6 [bdt = benzenedithiolato] complex, a synthetic mimic of the [FeFe] hydrogenase enzyme can electrochemically convert protons into molecular hydrogen. Molecular understanding of the cascade of reaction steps is important for the design of more eﬃcient catalysts. In this study, we investigate the reaction mechanism of the hydrogen production catalysis in explicit solution of acetonitrile using first principles molecular dynamics simulations. We have characterized all reduction and protonation intermediates taking part in the catalytic cycle. Free energy surfaces of the activated reaction steps are calculated using metadynamics. We find that the second protonation leading to molecular hydrogen formation is the rate limiting step. Direct protonation of the bridging hydride by a proton from the solution to form H 2 is the most favorable reaction pathway. However, also a bdt sulfur atom can become protonated, leading to a possible proton trap state that reduces the catalytic eﬃciency. Our calculations validate the ECEC mechanism proposed using cyclic voltammetry.


Introduction
Molecular hydrogen is a promising alternative source of sustainable energy 1 and the development of inexpensive catalysts for hydrogen fuel production is therefore a research topic of paramount importance. A particularly interesting route to developing efficient catalysts using earth abundant metals is inspired by enzymes that catalyse hydrogen formation. Hydrogenases form a class of enzymes found in nature, which can reduce protons to form molecular hydrogen. 2 Hydrogenases are classified into three main categories based on the metal content in the enzyme: 3 [NiFe] hydrogenase, [FeFe] hydrogenase and Fe hydrogenase. Diiron hydrogenase is considered to be the most active enzyme for proton reduction with a turnover frequency of the order of B10 000 s À1 . 4 The active part of the enzyme is called the H-cluster (see Fig. 1). The H-cluster consists of a Fe 2 S 2 unit linked to a Fe 4 S 4 cubane cluster via a cysteine. 3 The Fe 2 S 2 group contains CO and CN ligands and the two S atoms of Fe 2 S 2 unit are coordinated by an azadithiolato ligand.
A large variety of structural and functional mimics of [FeFe] hydrogenase have been synthesized, 5 of which the majority fall in one of the following three classes: those with an aromatic bridging ligand (e.g. benzenedithiolato (bdt)), [6][7][8][9][10][11] with an alkylic bridging ligand (e.g. ethanedithiolato (edt) and propanedithiolato (pdt)) [12][13][14] and those with a bridging azodithiolato (adt) ligand [15][16][17][18] (see Fig. 2). Catalysts with the aromatic bridgeheads are particularly interesting because they can undergo a reversible twoelectron reduction at an electrode potential of À1.3 V (vs. the Fc + /Fc 0 electrode). This is 0.37 V more than that of the aliphatic (pdt) bridgehead ligand, which undergoes an irreversible reduction at À1.67 V. The two-electron reduced intermediate has a pK a less than 23 and can be protonated using weak acids. Derivatives of complex 1 have been used on electrodes for the hydrogen evolution reaction. The reaction proceeds through a  bridging hydride intermediate of the diiron catalyst, which is different from that of the natural hydrogenase enzyme, in which the reaction intermediate is a terminal hydride. It is speculated that this difference in the reaction pathway could be the reason for the lower activity of the synthetic mimics. 19 More complex ligands (instead of CO) have been used in the synthesis of mimics over the last decade, which led to improvement of the catalytic activity, 20 however still none of them are as active as the natural enzymes. In addition, design of a diiron hydrogenase inspired catalyst that is soluble and stable in water solvent is particularly challenging and a topic of intense research. A detailed understanding of the mechanism of proton reduction is crucial for improving the performance of synthetic catalysts.
Computational methods have been used to study the mechanism of [FeFe] hydrogenase enzyme and mimics. 6,[21][22][23][24][25][26] This includes optimising the structures of possible intermediates taking part in the reaction cycle and calculating their reduction potentials and acidities. Most of these studies aimed at unraveling the catalytic mechanism of the catalyst or the enzyme are performed either at the density functional theory (DFT) level of theory augmented with the COSMO implicit solvation model, 26,27 or using classical molecular dynamics 25 or a hybrid quantum chemical/forcefield (i.e. QM/MM) approach. 28 While recent studies have provided interesting insights into the catalytic cycle of the diiron complex 1, 6,26 we still lack a complete understanding of the reaction mechanism and the redox properties of the intermediates. Fig. 3 depicts the proposed reaction mechanism for proton reduction, catalysed by compound 1. Cyclic voltammetry and computational studies confirm that 1 undergoes a two electron reduction to form 1 2À . The protonation of 1 2À will yield 1H À and further reduction of the 1H À intermediate results in hydrogen evolution in the presence of an acid with pK a o 23. 6 In our present work, we performed an extensive theoretical study using first principles molecular dynamics simulations of complex 1 explicitly solvated in acetonitrile. The geometries and electronic structures of all the different intermediates in the proposed reaction mechanism are discussed. Energy barriers for protonation and hydrogen formation are calculated using the metadynamics enhanced sampling method. 29,30 Prediction of formal (integer valued) oxidation states from the quantum chemical calculations is not straightforward. 31 Population analyses [32][33][34] can give useful information about the charge distribution in the complex, however the computed partial charges generally do not provide a clear insight into the formal oxidation state of metal centers or larger complexes. Instead here, we determine the oxidation states using so-called maximally localized Wannier functions. 35,36 This approach provides a clear understanding of the reaction mechanism, including the free energy barriers of the important reaction steps and the oxidation states of the intermediates.

Methods
All DFT-MD simulations were carried out using the mixed Gaussian and plane wave method as implemented in the CP2K package. 37 We used the PBE exchange correlation functional 38 augmented with Grimme's D3 39 dispersion correction and the DZVP-MOLOPT 40 Gaussian basis set with a plane wave cut off of 300 Ry. GTH-type pseudo potentials 41,42 were used to represent the valence-core interactions.
The molecular system contained complex 1 solvated in 40 acetonitrile molecules in a cubic supercell with a length of 16.2163 Å subject to periodic boundary conditions. The density was determined from a classical molecular dynamics simulation in the NPT ensemble. All following DFT-MD simulations were performed in the NVT ensemble with a CSVR 43 thermostat with a time constant of 50 fs maintaining an average temperature of 298.15 K. All DFT-MD simulations of reaction intermediates of complex 1 were carried out in the electronic low-spin state, which we established to be the ground-state with DFT geometry optimisation calculations (see Fig. 1 and Tables S1-S6 in the ESI † for details). A DFT optimised neutral structure was solvated in acetonitrile solvent and equilibrated with DFT-MD. After the equilibration run, electrons and protons were subsequently added to the system to create the relevant catalysis intermediates. DFT-MD simulations were carried out for 7 ps for each of the intermediates described in the Results section, of which the first 2 ps were used for equilibration.
Wannier center analysis 44 as implemented in CP2K was used to track the changes in the oxidation states. An electron counting technique was used to assign the oxidation state, as follows. [45][46][47] For a closed shell calculation, each Wannier center corresponds to two electrons. If the Wannier center is either inside the atom or very close to one of the atoms in a covalent bond, we assign the Wannier center to that atom. If the Wannier center is at the middle of a covalent bond, we assign one electron to each of the atoms. For an open shell system, the same scheme is used but with each Wannier center corresponding to only one electron.
Metadynamics as implemented in PLUMED 48 was used for the free energy calculations of the protonation reactions. We used   Gaussian potentials with a height of 0.5 kJ mol À1 and widths of 0.1 Å, with a deposit interval of 50 fs. Avogadro 49 was used for setting up the initial structures. The VMD 50 tool was used for visualisation.

Results and discussions
Below, we will discuss first the characterisation of the different intermediates taking part in the reaction mechanism shown in Fig. 3, including the calculation of the formal oxidation states.
In the following section, we present the first protonation reaction step and the calculation of the pK a . The final section presents the energy barriers for the hydrogen formation.

Characterisation of compound 1 in different oxidation and protonation states in acetonitrile solution
Each of the six reactant, intermediate, and product species occurring in the catalytic cycle starting from compound 1, shown in Fig. 3, was simulated and characterized in an equilibrium NVT DFT-MD simulation at T = 298.15 K. Representative snapshots of the diiron complexes from these simulations are shown in Fig. 4 for the unprotonated species and in Fig. 5 for the protonated complexes.
Starting from the neutral complex in Fig. 4a, the complex maintains its highly symmetric structure, known from static DFT geometry optimisations, throughout the simulation. The complex appears particularly stable, with none of the four Fe-S bonds showing elongated vibrations and all CO groups remaining in terminal positions, rather than a bridging configuration. Average bond lengths and angles are in agreement with crystallographic data. We use the Wannier center analysis to obtain formal charges and oxidation states of the ligands and iron centers. In Fig. 4, the Wannier centers that contribute to the electron count of the iron ions are shown as orange balls. In the neutral complex, seven Wannier centers are located at each iron and one Wannier center is seen in the middle between the irons. For this closed shell calculation, each Wannier center represents two electrons, and the iron core pseudo-potential has a positive charge of 16, so that the formal oxidation state of each iron is +1. The two sulfur atoms each have a formal charge of À1 and the carbonyl ligands are neutral. This Wannier center analysis assignment is in agreement with experimental results. [51][52][53][54] Reduction of 1 is a two electron reduction process and the 1 À intermediate reduces instantaneously further to 1 2À , which makes it difficult to characterise the singly reduced species experimentally. One electron reduction of complex 1 changes the geometry of the complex. A detailed study of this intermediate 11 shows all the CO ligands in terminal positions and an elongated Fe-S bond. A previous computational study 6 suggested a bridging CO ligand for this intermediate. The energy difference between the two isomers is about 3 kJ mol À1 . 11 Fig. 4b shows the structure of the 1 À intermediate in acetonitrile solution taken from our DFT-MD simulation. We observe that Fe-S bond cleavage occurs instantaneously after the first one electron reduction. All CO groups remain in terminal positions during the simulation, in agreement with experimental observations. 11 Since this is now an open shell calculation (the spin multiplicity of the complex is 2), one Wannier center represents one electron. The number of Wannier centers for the two Fe atoms are 8 and 7 respectively in the spin up configuration, and 7 and 7 for the spin down configuration. There is one Wannier center in the middle of the two Fe atoms for each spin; we assign one electron to each Fe center. This yields a total of 16 and 15 electrons respectively and the calculated oxidation states are thus 0 and +1. This again matches very well with experiments. [51][52][53] The second reduction gives the 1 2À intermediate, which is seen to deform even further from the symmetric initial neutral complex as seen in Fig. 4c. In the simulation, Fe-S bond breaking and simultaneous rearrangement of a terminal carbonyl group to the bridging position takes place instantaneously. After this initial reconstruction, the structure remains stable throughout the remainder of the simulation. The Wannier center analysis shows eight centers (representing 2 electrons each) at each Fe atom. The calculated oxidation state of the two Fe atoms is therefore 0. The 1H complex is formed when the 1 À intermediate is protonated by the presence of an acid with a pK a r 12.7. Mirmohades et al. 11 characterized 1H experimentally and our simulation results are in agreement with their study. 1H has a symmetric structure similar to 1, where the Fe-S bonds are intact and all the CO ligands are at terminal positions. The hydride is at a bridging position in between the two iron atoms, as shown in Fig. 5a. The oxidation state is calculated as HFe II -Fe I . There is one Wannier center located at the H atom at each spin state making a total of two electrons at the hydrogen, which yields a formal charge of 1 À for the hydrogen, a clear confirmation of the bridging hydride species. Second protonation of this intermediate is only favorable at a lower pH and the pathway toward molecular hydrogen formation is through further reduction of this intermediate to the more active HFe I -Fe I complex. 53 For the protonated doubly reduced complex, 1H À , we characterise the two most favorable isomers: a bridging hydride intermediate and a sulfur protonated intermediate. The geometry of the 1H À bridging hydride intermediate is shown in Fig. 5b. Note that it has both a bridging carbonyl ligand and a bridging hydride. The formal charge of each the two Fe ions is +1 according to our Wannier center analysis. The HFe I -Fe I intermediate is postulated to be an important intermediate in the catalytic cycle. [51][52][53] The doubly reduced complex may also be protonated at one of the sulfur atoms, instead of at an iron. This 1SH À species (Fe 0 -Fe 0 ) forms a hydrogen bond with a solvent acetonitrile nitrogen atom during the simulation, see Fig. 6 for an illustrative snapshot. In case of the bridging hydride isomer, no hydrogen bond formation between the hydride and the solvent was observed. Although less favorable than the bridging hydride isomer, no tautomerization via proton transfer was observed during the DFT-MD simulations.
Cyclic voltammetry studies of 1 have shown that one electron reduction of 1H À takes place to form 1H 2À and this reduction leads to formation of molecular hydrogen. 6 The structure of 1H 2À shown in Fig. 5b is similar to that of 1H À . During the initial stage of the DFT-MD simulation after reduction of 1H À , we observe another CO group moving to the bridging position. There are thus temporary two bridging CO groups and a bridging hydride and no bridging sulfur atoms. But this structure is not very stable and rearranges to a single bridging CO and a bridging sulfur. A previous computational study also reported a bridging ligand and a bridging hydride 6 Fig. 7 provides an overview of the iron oxidation states in each of the intermediates taking part of the catalytic cycle.

Free energies of the protonation reactions
The protonation reaction free energies of reaction steps 1 À + H + -1H and 1 2À + H + -1H À (see also Fig. 3) are directly related to the pK a values of the protonated species through eqn (1): in which DG is the reaction free energy, R is the ideal gas constant, and T is the absolute temperature. We use a series constrained This journal is © the Owner Societies 2020 Phys. Chem. Chem. Phys., 2020, 22, 10447--10454 | 10451 molecular dynamics simulations, following the example of ref. 55, to compute the free energy profile as a function of a simple geometric reaction coordinate. Rather than a holonomic constraint, we use a stiff harmonic spring, with a force constant of k = 100 kJ mol À1 Å À2 to restrain the sampling at specific values of the reaction coordinate. Although under experimental conditions the proton is provided by a suitable acid molecule, here we use a solvent acetonitrile molecule as the proton donor to directly obtain pK a values corresponding to this solvent. 3.2.1 pK a calculation of 1 À . We have calculated the free energy profile for the protonation reaction of 1 À : 1 À + H + -1H, using a series of 9 constrained DFT-MD simulations. The reaction coordinate was the distance between the proton and the nitrogen of the donating acetonitrile molecule. An equilibration run of 3 ps and a production run of 5 ps were carried out at each constrained N-H distance value. Integration of the sampled average constraint force resulted in the free energy profile shown in Fig. 8 and an estimate of the protonation free energy of DG = À68 kJ mol À1 . Using eqn (1) to compute the pK a of 1 À results in a value of 12, which is close to the experimental value of 12.7. 11 Protonation of catalyst 1 in the singly reduced state thus requires a rather strong acid (note that pK a values are shifted in acetonitrile with respect to those in water solvent, for example for carboxylic acids the shift is on average about 15.5 pK a units).
3.2.2 Protonation of 1 2À . Although the previous protonation reaction of intermediate 1 À was found to be a downhill process with a minuscule free energy barrier, we did not observe a spontaneous protonation during an equilibrium DFT-MD simulation starting from a configuration with a nearby protonated acetonitrile molecule. Instead for the doubly reduced system, the protonation of 1 2À was seen to take place on the picosecond time scale of an equilibrium DFT-MD simulation, leading to the bridging hydride configuration, as illustrated by snapshots from the initial and final configurations in Fig. 10. Both experimental 56 and theoretical 57 studies have reported that the bridging position is the most stable position for the first proton. Our simulation confirms that protonation of the doubly reduced complex is a spontaneous low-barrier process and not a rate limiting step for the reaction mechanism. 8 Experiments have shown that this protonation requires a relatively weak acid with an acidity constant of pK a r 23. We have calculated the pK a of this reaction using constrained DFT-MD simulations. Fig. 9 shows the computed free energy profile. The resulting reaction free energy is 116 kJ mol À1 , which corresponds to a pK a of 20.3, in good agreement with the experimental observation.
In Section 3.1, we noted that the 1 2À intermediate can also be protonated at one of the sulfur atoms of the bdt ligand, forming the 1SH À complex (see Fig. 6), which is a tautomer of the bridging hydride isomer of 1H À . For completeness, we also performed a DFT-MD simulation starting from a configuration with a protonated acetonitrile molecule closer to a sulfur atom. An immediate transfer of the proton from the solvent molecule to the bdt sulfur occurs. Sulfur protonation has been a topic of research using computational methods, 21,57 but there is no experimental characterisation of this intermediate to the best of our knowledge. The bridging hydride is the most stable isomer. 56,57

Metadynamics simulation of the H 2 formation
The final reaction step in the catalytic cycle is the formation of molecular hydrogen by (second) protonation of the, now, 1H 2À intermediate. We investigated two possible pathways for this step: (1) direct protonation of the bridging hydride by proton donation from a donor in the solvent, and (2) a two-step process via a doubly protonated intermediate that has a protonated sulfur atom and a bridging hydride. Because the protonation   This journal is © the Owner Societies 2020 reaction and subsequent H 2 release involves forming and breaking of three or more bonds, computing the free energy profile along a single reaction coordinate with constrained DFT-MD cannot be done. Instead, we use metadynamics in combination with DFT-MD to compute the free energy landscapes as a function of two reaction coordinates, or collective variables.
The metadynamics simulation was preceded by a DFT-MD equilibration simulation starting from a protonated solvent acetonitrile molecule close to the complex' bridging hydride. No spontaneous proton transfer was observed, suggesting that this reaction is activated and a rare event on the picosecond time scale of our simulations. For the metadynamics simulation of the direct protonation mechanism, we employed two collective variables to describe the reaction: (1) the difference of the distance between the bridging hydride and the center of mass of the two iron atoms, d(C-H1), and the distance between the bridging hydride and the second proton, d(H1-H2); and (2) the difference of the distance between the solvent nitrogen atom and the (second) proton, d (N-H2), and the distance between the bridging hydride and the second proton, d (H1-H2). The metadynamics simulation is stopped after the reactant free energy well is estimated and the barrier is crossed, and the formed H 2 molecule diffuses into the solvent. The resulting free energy surface is shown in Fig. 11, and the free energy barrier is 20 kJ mol À1 .
The alternative two-step mechanism of the H 2 formation reaction proceeds by first forming a doubly protonated intermediate by protonation of a sulfur atom. The metadynamics result of this first step is shown in Fig. 12. The two collective variables were: (1) the distance between the sulfur atom and the (second) proton, d(S-H), and (2) the distance between the acetonitrile nitrogen atom and the proton, d(N-H). The protonation