Successes, challenges, and opportunities for quantum chemistry in understanding metalloenzymes for solar fuels research ChemComm FEATURE ARTICLE

Quantum chemical approaches today are a powerful tool to study the properties and reactivity of metalloenzymes. In the field of solar fuels research these involve predominantly photosystem II and hydrogenases, which catalyze water oxidation and hydrogen evolution, as well as related biomimetic and bio-inspired models. Theoretical methods are extensively used to better comprehend the nature of catalytic intermediates, establish important structure–function and structure–property correlations, elucidate functional principles, and uncover the catalytic activity of these complex systems by unravelling the key steps of their reaction mechanism. Computations in the field of water oxidation and hydrogen evolution are used as predictive tools to elucidate structures, explain and synthesize complex experimental observations from advanced spectroscopic techniques, rationalize reactivity on the basis of atomistic models and electronic structure, and guide the design of new synthetic targets. This feature article covers recent advances in the application of quantum chemical methods for understanding the nature of catalytic intermediates and the mechanism by which photosystem II and hydrogenases achieve their function, and points at essential questions that remain only partly answered and at challenges that will have to be met by future advances and applications of quantum and computational chemistry.


Introduction
The pressing need to develop energy sources and alternative fuels that are based on renewable, environmentally friendly, and affordable approaches is critical for successfully facing the mounting energy and climate challenges. 1,2 Molecular solar fuels, such as hydrogen, 3 that can in principle be obtained by mimicking biological processes, represent the most promising a Aix-Marseille Université, CNRS, iSm2, Marseille, France. E-mail: maylis.orio@univ-amu.fr b Max-Planck-Institut für Kohlenforschung, Kaiser-Wilhelm-Platz 1, 45470 Mülheim an der Ruhr, Germany. E-mail: dimitrios.pantazis@kofo.mpg.de

Maylis Orio
Maylis Orio was born in 1981. She obtained her PhD in Chemistry from the University of Grenoble in 2007. She is a CNRS researcher from Aix-Marseille University. Her research interests focus on the use of quantum bioinorganic chemistry in relation with catalysis, including in the field of hydrogen production.

Dimitrios A. Pantazis
Dimitrios A. Pantazis was born in 1976. He obtained his PhD in Chemistry from the University of York and is currently a research group leader at the Max-Planck-Institut für Kohlenforschung. His research interests focus on quantum bioinorganic chemistry and theoretical spectroscopy, particularly in relation to natural and artificial photosynthesis.
way of meeting this need. In the hydrogen evolution reaction (HER), electrons and protons are combined to generate molecular H 2 , which stores energy in its chemical bond. Hydrogen production through water splitting is one of the preferred solutions in the long run for the storage of renewable energy. Today hydrogen is principally used for industrial processes. However, it has the potential to play a major role as energy carrier in transportation, the gas industry, as well as in the generation of electricity and heat. In addition to other technical challenges, for the expansion of H 2 usage in these directions it would be necessary to increase the global production of H 2 considerably and sustainably. 4 Nature can provide useful guidelines in this quest, because metalloenzymes involved in water oxidation and hydrogen evolution serve as archetypes for these reactions and for the development of artificial catalysts for solar fuels. 2,5-14 Nature offers a unique example of a water-oxidizing enzyme, photosystem II (PSII). 15,16 This harbours the oxygen-evolving complex (OEC) with its manganese-calcium cluster (Scheme 1a) that catalyses the oxidation of water into dioxygen, protons, and electrons, powered by light-driven charge separation. Several aspects of the function of PSII remain unclear, 17 while the development of molecular manganese-based structural mimics of the OEC [18][19][20][21][22][23][24][25] has not been accompanied by comparable catalytic water oxidizing activity. [26][27][28] In terms of H 2 evolution, biology provides good examples in hydrogenases (HG). These enzymes contain nickel and iron active sites (Scheme 1b) and achieve catalytic performances that rival synthetic platinum-based catalysts for hydrogen production. [29][30][31][32] Therefore, in analogy to the manganesebased water oxidizing PSII, hydrogenases serve to inspire the design of artificial hydrogen evolution catalysts based on Earthabundant metals. 33 Several studies have reported homo- [34][35][36][37] and hetero-polynuclear [38][39][40][41] metal complexes as electrocatalysts for HER reproducing the structure of the active site of hydrogenases. In addition to biomimetic approaches, a variety of molecular electrocatalysts have been developed to facilitate hydrogen production by making it occur faster with a lower energy input. Up to now, platinum-based catalysts provide the best performance in hydrogen evolution. 42 However, due to its scarcity and cost, intense research efforts have been devoted in developing non-noble transition metal catalysts. 43 In this connection, diverse ligand skeletons with cobalt, iron, and nickel complexes [44][45][46][47][48][49][50][51] were shown to evolve hydrogen both in organic and aqueous media (Scheme 2).
Research in the PSII and HG fields has benefited greatly from the contributions of quantum chemistry. During the last two decades of the 20th century the density functional theory (DFT) revolution extended the applicability of quantum chemistry to realistic transition metal systems in general and metalloenzymes in particular. [52][53][54][55][56] Particularly after the advent of hybrid functionals it was possible to predict molecular structures and reaction profiles with sufficient accuracy to make meaningful analysis of geometric and mechanistic aspects of bioinorganic systems. The further combination of DFT with classical force-field methods within QM/MM approaches 57 has become a considerably powerful component of modern computational chemistry.
The first two decades of the 21st century broadened even more the realm of applicability of quantum chemistry to metalloenzymes. One major development is that ''quantum chemistry'' in this field is no longer synonymous with ''DFT'', but relentless advances in methodology, implementation, and infrastructure have made single and multireference wave function theory a practical and accessible option. Moreover, quantum chemistry today reaches much further than mere geometry optimizations and simple reaction energetics. Theoretical advances mean that much deeper insights into the electronic structure and properties can be obtained, and hence much deeper and more insightful connections with experiment can be established.
The two title systems of the present feature article, photosystem II and hydrogenases, have served as application platforms on which quantum chemistry has reached new heights in terms of what can be computed, and how it advances the broader research fields. In the following we will describe selected 2. Photosystem II Photosynthesis involves light harvesting, excitation energy transfer, and creation of charge separated states that drive subsequent redox transformations. 15,16 In oxygenic photosynthesis the pigment-protein complex PSII utilizes this charge separation, likely initiated at the primary donor Chl D1 , [64][65][66][67] to create the most oxidizing agent in biology, a radical cation (P680 + ) delocalized among a pair of chlorophylls (P D1 /P D2 ). This oxidizes the Mn 4 CaO x cluster of the OEC in successive one-electron abstraction steps via an intermediary redox-active tyrosine (Y Z ), until the OEC becomes able to complete the four-electron oxidation of water to dioxygen. The intermediate ''charging'' or ''storage'' states of the OEC are known as S i states, where the subscript denotes the number of stored oxidizing equivalents (i = 0-4). S 1 is the resting (darkstable) state with commonly accepted Mn oxidation states of Mn(III) 2 Mn(IV) 2 , while S 4 is an unobserved transient state that evolves dioxygen (Fig. 1). Electron and proton transfers take place alternately, 68 maintaining a narrow potential range for the OEC for all S i -S i+1 state transitions. On the electron acceptor side of PSII a plastoquinone is doubly reduced and released from the complex to transfer reducing equivalents further along the photosynthetic chain, eventually to be used in CO 2 fixation.
Our primary focus here is on the site of water oxidation and the contributions of quantum chemistry in understanding its nature, its properties, and its function. We will cover selected developments and applications in recent years, mostly from our group, that had significant impact on discovering specific correlations between structural features and spectroscopic properties of the OEC, and which contributed to interpreting experimental results and directing further investigations. We will also identify outstanding challenges that will require more than simple extensions of present quantum chemical approaches to be successfully met in the future.

Correlating molecular structure with spin states
For a long time geometric information on the OEC was confined to metal-metal distances derived from EXAFS. 69 This allowed discussion of plausible structural features but could not lead to well-defined atomistic model that are required for quantum chemical investigations. Even the first crystallographic models [70][71][72] were not of sufficient resolution to allow definition of a unique structure. Therefore, early computational studies 73,74 explored a range of geometric and mechanistic possibilities but incorporated numerous structural assumptions that later proved inaccurate. Crystallographic models of sufficiently high resolution that appeared in the last decade and that more recently utilize pulses from X-ray free electron lasers changed this situation by accurately defining the spatial arrangement of metal ions and oxygen bridges as well as the type and mode of coordination of first and second sphere ligands. [75][76][77][78][79][80][81] As shown in Scheme 1a the core of the OEC consists of four Mn and one Ca ions connected by oxo bridges and additionally bridged mostly by carboxylate residues, while there are also two terminal water-derived ligands (H 2 O or OH) on each of Mn4 and Ca. The redox-active tyrosine residue (Tyr161, or Y Z ) is in close proximity and interacts with Ca and the rest of the inorganic cluster via hydrogen bonding. Ambiguities remained, 82 and they still do, 17 for example regarding the protonation states of oxo bridges and the extent of radiation-induced reduction of metal ions, 83-85 but these advances nevertheless enabled quantum chemical studies to be performed on a much safer basis than ever before. 86 One important consequence was that the spatial arrangement was sufficiently well defined to allow theoretical studies of electronic structure, spin states, and spectroscopic properties beyond superficial energy comparisons, thus enabling for the first time direct connections between quantum chemical studies and the vast repository of spectroscopic observations on various catalytic intermediates.
The most obvious target property to begin with is the spin state of specific catalytic intermediates, which in many cases had been identified, or at least constrained, by EPR spectroscopy. Advances in this direction were tied to the investigation of magnetostructural correlations and began with the spectroscopically well-studied S 2 state of the OEC and synthetic analogues. The S 2 state was long known to exhibit two types of EPR signal that correspond either to a low-spin (S = 1/2, g E 2) multiline state or to a higher-spin state (S Z 5/2, g Z 4.1). [87][88][89][90][91][92] A landmark study in 2012 established a correspondence between two structural forms of the OEC core and these two spin states. 93 The structural forms were derived from refinement of the 2011 crystallographic model. 75 They were almost isoenergetic and featured different distribution of Mn valences (III-IV-IV-IV versus IV-IV-IV-III) and slightly different connectivity of an oxo bridge that led to their description as open-cubane (S 2 A ) and closed-cubane (S 2 B ). 93 As will be explained below, quantum chemical investigations of the spindependent properties of these models supported their association with the two observed EPR signals.
The quantum chemical approach used for determining the spin states is based on a generalization of the two-centre broken symmetry (BS) DFT methodology [94][95][96][97] to the tetranuclear case. The low-lying states of a cluster representing the OEC can be described using the effective Heisenberg-Dirac-van Vleck Hamiltonian (HDvV), which describes magnetic levels in terms of fictitious interactions between local spin magnetic moments: where Ŝ i and Ŝ j are the spin operators and J ij are the pairwise exchange coupling constants that parameterize the fictitious spin interaction. The sign of J ij signifies whether the interaction is antiferromagnetic ( J ij o 0) or ferromagnetic ( J ij 4 0). Diagonalization of the Hamiltonian yields the spectrum of spin states (''spin ladder'') that can be directly compared to experiment as long as the effective Hamiltonian is a valid representation of lowenergy states. The broken symmetry approach attempts to determine the J ij parameters employing a set of Kohn-Sham determinants in which the unpaired electrons are localised at the Mn ions in all possible up/down (a/b) combinations while maintaining the local high-spin electronic configuration on each site. Using as an example the S 2 state of the OEC with Mn oxidation states III-IV-IV-IV, local S z of 2 for the d 4 ion Mn(III) and 3/2 for the d 3 ion Mn(IV), the high-spin determinant can be described as |+2, +3/2, +3/2, +3/2i (or ''aaaa'', M S = 13/2) and BS determinants are created by inverting (''flipping'') local spins of Mn ions. For example, |+2, +3/2, +3/2, À3/2i (''aaab'', M S = 7/2) or |+2, À3/2, À3/2, +3/2i (''abba'', M S = 1/2). Each BS determinant is characterized by an M S value and not S, because it is not a spin eigenfunction. The energies of the BS determinants can be used to deduce a set of J ij values for the S z -only (Ising) Hamiltonian. Subsequently, the same exchange coupling constants are used to diagonalize the full HDvV Hamiltonian and obtain the spin ladder. The details of the procedure have been described extensively in the literature. [98][99][100] Fig. 2 provides a schematic flow diagram that summarizes the points made above. This lengthy procedure obviously involves significantly increased computational and human effort compared to studies that simply ignore the spin state problem. For each structural model to be considered one has to complete a series of calculations that often exhibit challenging convergence behaviour, and then perform additional analysis to extract exchange coupling constants and magnetic spin states. Yet, despite the considerable overhead, the results have been richly rewarding and have pushed research in natural water oxidation in new and fruitful directions. It is crucial to stress that no discussion of spin states can be made by comparing energies of different broken-symmetry determinant, simply because these do not correspond to observable states. It is an unfortunate fact that several papers erroneously discuss spin states by direct reference to broken-symmetry solutions.
The application of the approach to the two valence isomers found for the S 2 state of the OEC showed that they differ in their exchange coupling scheme (Fig. 3). 93 The isomer labelled as S 2 A (III-IV-IV-IV) has antiferromagnetic interaction between Mn1-Mn2 and Mn3-Mn4 resulting in a spin S = 1/2 ground state, Fig. 2 Flow diagram of calculations in the context of the brokensymmetry DFT approach for determining the exchange coupling constants and spin state ladder for a computational model of the S 2 state of the OEC, assuming Mn oxidation states of III-IV-IV-IV. Note that in this case there are only seven BS determinants that can be formed in addition to the highspin determinant, yet the proper solution to the exchange coupling problem yields in total 320 states (including sublevels of spin multiplets). whereas S 2 B features strong ferromagnetic interactions within the cuboidal Mn1-Mn2-Mn3 unit (effective spin of the cube S = 9/2) resulting in a total ground state spin of S = 5/2 for this isomer. This suggests already a potential correspondence between the two structural models and the observed spin states. This was further supported by explicit simulations of EPR spectra, which associated the S 2 A with the g = 2 EPR signal and the S 2 B model with the g = 4.1 signal. 93 Importantly, the results on the magnetic properties of the S 2 state isomers are aligned with experimental and theoretical results on synthetic trinuclear and tetranuclear analogues. 101,102 For example, synthetic Mn(IV) 3 Ca cubane complexes 22,103 that mimic this structural unit of the S 2 B model similarly have dominant ferromagnetic coupling among their Mn ions and display high-spin (S = 9/2) ground states. 102 This can be explained at the orbital level by elimination of superexchange interactions between the Mn ions due to the angles of the oxo bridges. 102 The broken-symmetry approach can be equally well applied to spinfrustrated systems with antiferromagnetic ground states, as exemplified by a combined experimental and theoretical study of a lowspin trinuclear terpyridine cluster with a Mn(IV) 3 O 4 core. 101 Consideration of spin state alone may be insufficient to establish a firm one-to-one correspondence between a computational model and a true catalytic intermediate. For example, other structural models or protonation patterns can be shown to have the same spin state as that observed experimentally. 104,105 It is only the combined evaluation of geometries, spin states, energetics, and spectroscopic parameters that can lead to safe assignments. In any case, spin states are strong constraints and essential for discussing the relevance of any computational model to an experimentally observable state of the physical system. Thus, they have been central in many studies that examined the nature of specific intermediates, by showing connections between structural features and magnetic properties. For example, BS-DFT analysis of spin states for various isomeric forms of the S 3 state established that the S = 3 states observed by magnetic resonance spectroscopies 106-108 must be attributed to structural forms that have one additional water-derived ligand compared to the preceding S 2 state, and hence forms where all four Mn ions are sixcoordinated. 104,106 By contrast, a form of the S 3 state that has no additional ligand and contains a five-coordinated Mn(IV) ion is characterized by the highest possible spin for a Mn(IV) 4 cluster, of S = 6 ( Fig. 4). 108,109 Similar structure-magnetism correlations have been used to examine the effect of protonation of oxo bridges on specific exchange coupling pathways. It is known that in general an oxo bridge enables stronger superexchange and thus enhances antiferromagnetic coupling between two Mn ions compared to a hydroxo bridge. [110][111][112] Applied to the OEC, study of the effects of protonation or hydrogen bonding to specific oxygen bridges contributed, among others, to screening of structural models for the S 0 state 104,113 and to evaluation of possibilities for the interaction of small molecules such as methanol and ammonia with the manganese cluster. [114][115][116][117] It is also important to note that since different oxidation states of the Mn ions lead to different exchange interactions, analysis of magnetic properties in terms of exchange couplings and spin ladders is essential for evaluating redox-isomeric forms of specific S-state intermediates 104 and for evaluating competing oxidation state paradigms. 104,118,119  The valence isomers rationalize the two EPR signals associated with the S 2 state. In addition, they are very close in energy and interconvertible over a low barrier, as required experimentally, and they are each consistent with independent spectroscopic measurements for each g signal, including the 55 Mn hyperfine coupling constants and the distinct localization of the Mn(III) ion in each case. These successful applications should not obscure the fact that the spin-projected broken-symmetry DFT methodology, even if applied correctly, depends considerably on the nature of the functional, on the nature of the system under study, and on the validity of the effective spin Hamiltonian assumed to describe the system. High-valent/high-spin manganese systems are among the most favourable systems for this approach. Indeed, according to a study by Krewald et al., 104 the brokensymmetry approach using the TPSSh functional and the complete methodology for calculation of exchange coupling constants and spin states predicts with absolute accuracy the correct ground spin state of all known oligonuclear manganese complexes with oxidation states relevant to the OEC. Exceptions may be associated with limitations of the isotropic exchange model rather than with failures at the DFT level. Of course, even if the energetics obtained by the above procedure are useful, it is always desirable to have direct access to the spin states themselves, which can only be achieved by using multireference methods. For multi-electron systems such as those containing high-valent Mn ions, multireference approaches have found limited application so far. 120,121 Notable studies that are directly relevant for the OEC include the use of density matrix renormalization group (DMRG) based complete active space selfconsistent field (CASSCF) and N-electron valence state perturbation theory (NEVPT2) calculations on dinuclear Mn and Cr complexes. 112,122,123 Such calculations are far from trivial. They require careful study of convergence with respect to the active space and other technical parameters, and they present a series of requirements that must be satisfied for the results to be meaningful and reliable. As defined by Roemelt et al., 122 these requirements include the necessity for state-averaged orbital optimization over all spin multiplets that are part of the exchange coupling problem and the necessity to include dynamic correlation. Results obtained without state averaging, without complete orbital optimization, and without additional treatment of dynamic correlation are numerically inadequate and may show substantial deviations from the correct order and spacing of spin states. In cases where such deviations have been reported for non-manganese systems, 124,125 it is not always clear whether the correct application of the methodology as defined above would suppress or eliminate them. 126 Systems with more than two Mn ions present profoundly greater difficulties for such methods. This is not strictly because of the larger active space but because of how the calculations should be performed in order to enable us to derive and discuss magnetic properties. Although state-specific DMRG-based CASSCF calculations have been demonstrated for tetramanganese models, 127,128 these calculations have not dealt with, and cannot address, the questions of exchange coupling and spin state energetics; it cannot even be verified in principle that the correct ground state is computed. The above stated requirements derived from application of DMRG-CASSCF/NEVPT2 on a dinuclear Mn complex, 122 particularly the need for state-averaged orbital optimization over all magnetic spin levels, imply a steep increase of computational complexity and cost for exchange coupling analysis of a trimanganese system, which is an open challenge for multireference treatments, and probably render the tetranuclear case intractable with standard algorithms on conventional computers.

Spectroscopic parameters
In addition to spin states, an important point of contact between quantum chemistry and experiment are the various forms of spectroscopy used to study catalytic intermediates of the OEC. Here we focus only on parameters relevant to magnetic resonance spectroscopies. Spectroscopy-oriented quantum chemistry is more than a powerful tool for analysing experimental results and more than a way of maximizing the extraction of information from experimental spectroscopic data. It allows access to information that lies beyond the limitations of specific experimental approaches and can achieve unification of disparate sources of information into a common interpretation that connects experimental observations with the electronic structure of computational structural models.
Hyperfine coupling constants of Mn ions contain rich information on the local and global electronic structure of the OEC. Here the great challenge for quantum chemistry has been how to properly compute Mn HFCs for an oligonuclear exchange coupled system with DFT. This was addressed with the development of a generally applicable theoretical methodology that allows the use of the standard broken-symmetry DFT approach and makes connections between the various levels of Hamiltonian approximations involved in the fitting of HFC parameters from experimental spectra and in the calculation of HFCs from approximate DFT. 98 The interested reader is referred to the original publication for methodological details. 98 It is noted that for comparisons to experimental values a scaling factor has to be used for the final spin projected HFCs because of the systematic underestimation of core spin polarization by DFT. 129,130 This factor is not universal and may differ depending on details of the methodology such as the type of functional and the treatment of relativity, as well as on the reference set of compounds used to derive the required scaling. 98,99,110,[131][132][133] New theoretical developments in this area will be important for making the calculations of HFCs by DFT more robust. 134,135 Recent advances in the calculation of HFCs by correlated wave function methods such as coupled cluster theory 136 are also very promising but their applicability and reliability remain incompletely explored. 137 Following a pilot application 98 on a synthetic tetramanganese complex that is a Mn(III)Mn(IV) 3 analogue 138 of the S 2 state, broken-symmetry DFT was used in combination with the above method to compute spin states and hyperfine coupling constants of a wide variety of structural models that were considered at the time as geometric possibilities for the OEC. 99 This study showed that many ideas derived from EXAFS studies 139 as well as practically all models reported in computational studies until that time [140][141][142][143] were not consistent with the spectroscopic data on the S 2 state. These methods were essential in establishing that model S 2 A described in Fig. 3

corresponds very well with the 55 Mn
HFCs determined experimentally by ENDOR spectroscopy, [144][145][146][147][148][149][150] with one large, two medium, and one smaller hyperfine coupling constant. 93 These results demonstrated the importance of spectroscopy-oriented quantum chemistry in discriminating between structural possibilities and encouraged numerous subsequent applications to the OEC and other exchange-coupled oligonuclear systems, including refinements and extensions of the method. 131,[151][152][153][154] An obvious use of quantum chemistry in the study of 55 Mn HFCs is in the evaluation of different oxidation state possibilities, both in terms of absolute oxidation levels for the whole cluster and in terms of internal distribution of local oxidation states (valence isomeric forms) for a given total oxidation level or S-state intermediate. This approach has been adopted for evaluating various isomers of the half-integer-spin S 0 and S 2 states that have been studied extensively by ENDOR spectroscopy. [144][145][146][147][148] In particular for the S 0 state the situation arises where in principle the total number of unpaired electrons can be distributed among the four Mn ions with a formal Mn(III) 3 Mn(IV) or Mn(II)Mn(III)Mn(IV) 2 configuration. In combination with analysis of exchange coupling and spin states, computed 55 Mn HFCs have helped in evaluating various possible structure models. 104,113 Another important use of computed 55 Mn HFCs has been in the evaluation of possible models for the S 3 state. ENDOR-detected NMR studies on one of the components of this state showed that all Mn ions were electronic similar and isotropic. 106 Quantum chemical calculations helped to show that all-Mn(IV) cluster models where all Mn ions are octahedrally coordinated due to coordination of an additional water-derived ligand compared to the S 2 state are most consistent with the data. 104,106 This carries significant weight for the persistent question of whether the S 2 -S 3 transition involves metal-centred or ligand-centred oxidation. [155][156][157] The calculations firmly support the Mn-centred option because they demonstrate that the experimentally determined 55 Mn HFCs for the S 3 state are not reproduced with a Mn(III)Mn(IV) 3 cluster such as in the preceding S 2 state. Owing in part to its heterogeneous nature, as clearly revealed by various magnetic resonance spectroscopic studies 106-108,158-160 -though not resolved by currently available structural methods which can only provide spatially or temporally averaged representations -the composition of the S 3 state remains under investigation. 161-164 55 Mn HFCs have also been used in evaluating possible protonation states of oxo bridges and Mn-bound waterderived terminal ligands, and hence determine the most likely overall protonation level and pattern of the active site. 100,104 This analysis eliminated the possibility of bridge protonation in the S 2 state and led to the conclusion that the terminal Mn in the S 2 state is most likely ligated by one water and one hydroxo ligand. 100 The focus on spectroscopic properties enables direct screening of computational models against experimental constraints and therefore guides the study of the catalytic mechanism in ways that are otherwise impossible. Although emphasis in the above was placed on metal ions, ligand HFCs provide significant additional information. It is possible to use ligand HFCs (e.g. of Mn-coordinating oxygen or nitrogen nuclei) to assign metal oxidation states, 165 to study the identity and kinetics of exchangeable atoms, and to discover the mechanisms by which small molecules access the active site and interact with the Mn ions. 114,115,[166][167][168][169][170][171] An example that will be highlighted here concerns the interaction of the substrate analogue methanol with the OEC cluster. The importance of these studies is that they point to likely pathways for substrate delivery to the OEC. It has long been known that methanol modifies the magnetic energy levels and EPR signals of the manganese cluster of the OEC, but the access pathway(s) and mode of interaction remained unknown. 167,[172][173][174][175][176][177][178] Based on the experimental determination of the 13 C hyperfine parameters of isotopically labelled methanol interacting with the S 2 state of the OEC, 167 an extensive computational screening of structural models was carried out that involved calculation of spin state energetics and 13 C isotropic and dipolar HFCs for each one. 116 This led to rejection of several possibilities such as direct binding of methanol to Ca 2+ or one of the Mn ions, supporting instead a second-sphere interaction that resulted in reorganization of the hydrogen bonding network affecting the O4 bridge and the Mn3-Mn4 exchange coupling interaction (Fig. 5). 116 This is in contrast to the case of ammonia, which binds directly to Mn4, 115,170,171,179 at least in one of its interaction modes. The identification of a crucial difference between cyanobacteria and higher plants in a residue close to the proposed site of methanol binding (D1-N87A) 117 is consistent with the different response of plant versus cyanobacterial OEC to methanol, which can be explained by the restricted accessibility in the case of cyanobacterial OEC due to the more constrained channel architecture. 117 An important methodological advance in terms of studying and understanding the local anisotropy of the Mn ions was been the development of the local complete active space configuration interaction (L-CASCI) approach. 180 This enables the use of multireference methods for the calculation of local zero field splitting tensors and consists in conducting a multireference (CASCI) calculation with an active space specifically constructed to contain local orbitals for any given metal ion. Successive application of the method for each centre of interest eventually yields all on-site parameters for further analysis. Application of the L-CASCI approach to the S 2 state models of Fig. 3 allowed the explicit calculation of site ZFS tensors (specifically the spin-orbit coupling contribution D SOC ) for the Mn(III) ions in each valence isomer (Fig. 6). 180 The orientation of the D SOC principal axes correlate directly with the orientation of the pseudo-Jahn-Teller axes of the Mn(III) ions, which are oriented along the formally open coordination site of Mn(III) in each isomer. Importantly, the valence distribution, geometric parameters, and computed local ZFS for the Mn4(III) ion in the valence isomer S 2 B agree perfectly with analogous experimental studies of the g E 4.1 component of the OEC, 181 With D z oriented along the W2-Mn4-O5 axis. This strongly supports the assignment of S 2 B to the high-spin g E 4.1 signal, in contrast to other hypotheses that do not place the Mn(III) ion at the Mn4 site. The L-CASCI approach will certainly see increased use in the future as more in-depth analysis of catalytic intermediates is undertaken, not only for other intermediates of the OEC but also for other biological or synthetic clusters with multiple open-shell transition metal ions. Quantum chemical studies that focus on spectroscopic properties have also helped to elucidate critical aspects of the redox-active tyrosine Y Z and Y D residues. 182,183 Here it is the g-tensor of the tyrosyl radical that is the central spectroscopic parameter. This can be calculated with standard DFT approaches, although particular care has to be taken to ensure that the calculations correctly represent the electronic structure of both the Y Z radical and the Mn cluster. Most relevant for the OEC are studies that focus on the nature of tyrosyl radical intermediates formed during the S 2 -S 2 Y Z -S 3 transition. 109,182 Detailed investigations of the initial stages of the S 2 -S 3 transition showed that formation of the tyrosyl radical is accompanied by reorientation of the dipole moment of the OEC in the S 2 B state, 182 which likely triggers proton release from the terminal water ligand W1 of the Mn4 ion to the secondsphere Asp61 residue that acts as proton acceptor. [184][185][186] This deprotonation is necessary for subsequent oxidation of the Mn ion by the tyrosyl radical. 109 Y Z oxidation is accompanied by shifting of the phenolic proton to His190. The ultimate fate of this proton is still under investigation. [187][188][189] In terms of the spectroscopic properties of the formed tyrosyl radical, In the native system the g x value was relatively low compared to analogous literature precedents, which was attributed to the presence of three hydrogen-bonding interactions, two from adjacent water molecules and one from the His190 proton (Fig. 7). 182 Similar correlations derived from quantum chemical studies have been essential for probing the hydration and uncovering the structural basis of spectroscopic observations for the other redox-active tyrosine of PSII, Y D . 183 Intriguingly, the above observations suggest a direct role for Ca in structuring the hydrogen-bonding environment of Y Z and hence directly modulating its redox properties, 182 which is crucial for ensuring efficient coupling between the manganese cluster of the OEC and the charge-separation site (reaction centre) of PSII. Complementing this hypothesis, explicit calculations of the effect of calcium substitution by various redox-inactive cations showed that the redox potential of the OEC responds only to the total charge of the cation while remaining insensitive to its Lewis acidity, 190 in contrast to correlations established for synthetic complexes. [190][191][192]

Challenges in computational setup and modelling
The preceding sections highlighted contributions of quantum chemistry to understanding Nature's water oxidizing catalyst from the perspective of predicting magnetic and spectroscopic properties. However, it should be stressed that the use of an appropriate theoretical method or procedure to calculate geometries, energies, or other properties, in no way ensures attainment of reliable results that can be meaningfully correlated with experiment and used to evaluate computational models. This is because there are often practical aspects relating to the definition and treatment of the atomistic model itself that may introduce weaknesses which undermine the overall reliability of quantum chemical studies regardless of the ''intrinsic'' ability of any given theoretical method  to deliver high quality numbers. Unfortunately such weaknesses may remain unnoticed because the results obtained by calculations on inappropriate models are rarely, if ever, noticeably wrong in an absolute numerical sense. Here we single out a few issues in model setup that may compromise the relevance of computational results for understanding the real system. Although the examples are taken from OEC studies, the broader lessons are universally valid. When dealing with large and complex active sites such as the OEC, a principal concern is to build a quantum mechanics (QM) cluster model large enough to incorporate the effects of the immediate environment directly or indirectly. The QM model needs to be large enough to explicitly include the complete first coordination sphere and all elements of the second coordination sphere that interact via hydrogen bonding. Information about the structure of the protein matrix that is not explicitly included in the QM cluster model is included indirectly by constraining specific atoms. Constraints are typically applied to backbone atoms such as a-carbons and groups used to terminate peptide chains. When the QM cluster model omits hydrogens bonds that are present in the crystal structure, then additional geometric constraints may be considered, e.g. of specific dihedral angles, or at least optimized models should be examined very carefully for structural compliance with the absent protein matrix. A specific example of this is offered by His332, ligand to Mn1. Inspection of the complete environment of the OEC shows that the orientation of the His332 imidazole ring is restrained by a hydrogen bonding interaction with the peptide carbonyl of Glu329 (Fig. 8), which normally would not be included in QM cluster models of the OEC. The absence of this interaction may result in rotation of the His332 ring by ca. 901 during geometry optimization. Retegan et al. showed that such rotation directly affects the 55 Mn hyperfine parameters for Mn1 by modifying the nature of the Mn-N bonding. 193 Moreover, the incorrect rotation of His332 artificially increases the accessibility to Mn1, for example by water molecules. Therefore, the orientation of this residue affects considerably the computed relative energetics for water or hydroxy binding and the relevant models reported in the computational literature should always be scrutinized for structural correctness.
The combination of quantum mechanics with a force-field (molecular mechanics) treatment of the environment, i.e. the QM/MM approach, has a strong tradition in the OEC 140,186,[193][194][195] and it is expected to offer a better treatment of the environment when applied correctly. Here we would stress that the use of the QM/MM approach rarely allows for appreciable reduction in the size of the QM region compared to a standard QM cluster approach and does not negate the need to incorporate all important hydrogen bonding groups (including essential crystallographic waters) in the QM region. In fact, the use of a very small QM region leads with certainty to artefacts, which include spurious redox activity arising at the periphery of the QM region that results from redox imbalance (see below) or overpolarization from force-field charges.
Related to the above is the inclusion of the redox-active Tyr161 residue (Y Z ) in QM cluster models, or in the QM region of QM/MM models. Y Z is both structurally coupled to the rest of the OEC via a tight network of hydrogen bonds, and redoxcoupled as amply evident by the catalytic cycle depicted in Fig. 1 and by preceding discussions. Thus, omission of the Tyr161-His190 pair from computational models can lead to two critical artefacts: (a) removal of the hydrogen bonding network and drastic reorganization of the water molecules that are situated between Y Z and the Mn ions, and (b) absence of redox balance in the computational model, or inability to evaluate it. The first issue means that in the absence of the Tyr161-His190 pair and of the physiological associated hydrogen bonding network, calculations of hypothetical processes implicating these waters have unnaturally low barriers and little physical connection to the real system. The second issue means that without the tyrosine and without explicitly studying the relative energies of S i Y Z versus S i+1 Y Z states, no credible discussion of S-state transition energetics can be made. Explicit consideration of Y Z and of its redox role is not only essential when S-state transitions are examined or when Y Z participates directly in the mechanism. 109,182,[196][197][198] The redox properties of any computational model should be examined as standard practice to ensure that intrinsic redox balance is achieved and that no unphysical properties are inadvertently introduced. For example, a recent study showed that a set of assumptions regarding Mn oxidation states and protonation states of specific ligands results in the unphysical shifting of the site of subsequent oxidation from Y Z to the second-sphere histidine His337. 199

Catalytic progression and O-O bond formation
The precise details of the water oxidation mechanism during the later stages of the catalytic cycle remain a constantly evolving discussion. 17,[200][201][202] Open questions revolve around the steps of the S 2 -S 3 transition, the composition of the heterogeneous S 3 state and the ''active'' population of S 3 state that progresses further to form the O-O bond. Very little, and often indirect, experimental data exist on the S 3 Y Z intermediate and the most critical part of the cycle that involves dioxygen evolution and reconstitution of the S 0 state. [203][204][205] The two valence isomeric forms of the S 2 state of the OEC depicted in Fig. 3 were shown to behave differently upon   109,197 (see above for comments regarding the irrelevance of computational models that do not explicitly incorporate the Y Z -His190 pair and hence neglect the central question of redox balance between the manganese cluster and the redox-active tyrosine). Although various other possibilities may be envisioned, at this time the above S 2 -S 3 sequence appears to be the best attempt so far to accommodate most experimental observations on S 2 , S 2 Y Z and S 3 states, particularly from EPR spectroscopy.
In the above scenario a key species is the S 3 component with a five-coordinate Mn(IV) ion (S 3 B in Fig. 4), which must have a high spin state (S = 6) as predicted by the quantum chemical analysis of Retegan et al. 109 This has received very strong experimental support by a recent multifrequency EPR experimental study of the S 3 state in spinach PSII, which directly identified an S = 6 population with high effective anisotropy deriving from anisotropic exchange coupling between the coordinatively unsaturated Mn4(IV) ion and the Mn(IV) 3 Ca cubane subunit. 108 Importantly, the S = 6 component was found to be the exclusive species in methanol-treated spinach PSII and to form a major part of S 3 even in the native system. 108 Although the water bound intermediate-spin (S = 3) isomers of Fig. 4 have been experimentally confirmed as distinct populations of the S 3 state, Krewald et al. showed that water binding is not required to oxidize the cluster to the S 4 state. 206 Assuming the waterunbound form is catalytically active, this would allow Mn-based oxidation, leading to an S 4 state of the OEC with a five-coordinate Mn4(V) ion bearing a terminal oxo group at the formally ''W2'' position of earlier states. 206 This can then couple nucleophilically in an acid-base with the internal oxo bridge to form the O-O bond (Fig. 9a). 206 The above mechanism is the only currently discussed possibility that allows for complete separation of the charging phase (four Mn-cantered oxidations) from the actual catalytic phase, which can be achieved by even-electron O-O bond formation even in a single concerted step. 206 On the other hand, if a water-bound S 3 species is the active population of S 3 , then it is expected that the more stable S 3 A,W isomer will be involved. The reactivity in this case has been investigated by extensive computational studies and it is expected that the final S 3 -S 4 step involves ligand-cantered oxidation to create an oxyl radical. 143,207 This would subsequently couple in a radical oxyl-oxo coupling mechanism (Fig. 9b), which can itself be realized in multiple alternatives. 185,208,209 Note that both of the above scenarios involve homovalent Mn(IV) 4 S 3 isomers and assume a regular sequence of Mn-centered oxidations at least up to the S 3 state (for a radical oxyl-oxo coupling mechanism) or all the way to an S 4 state (for a nucleophilic oxo coupling mechanism). Several similar or distinct possibilities are being discussed, 81,161,162,164,[210][211][212][213] including alternatives avoiding alternative valence isomeric forms for the S 2 state, invoking early water binding before advancement to the S 3 state, assuming formation of an oxyl radical in the S 3 state or even ''early onset'' O-O bond formation already in the S 3 state, and others. Some of these alternatives may be in explicit or implicit disagreement with the most obvious interpretations of experimental data up to the S 3 state, 161 but in view of the absence of meaningful constraints past the S 3 state, they need to be carefully evaluated and for this reason further and more ambitious experimental and theoretical studies will be required.
A specific challenge for quantum chemistry in this case is to accurately distinguish energetically between different redox isomeric forms. Isobe et al. have already shown the sensitivity and precariousness of DFT when dealing with this problem. 163,211 The question of accurate energetics in this context is closely related to the general problem of spin-state energetics for transition metal systems, 214 though magnified in complexity by the presence of the four interacting open-shell Mn ions and the potential for oxygen-based radicals and O-O bonded intermediates. Achieving systematically accurate energetics in this case might be possible by a new generation of linear scaling coupled cluster approaches. A prominent example is the domain-based local pair natural orbital approach to the popular coupled cluster theory with single, double, and perturbative triple excitations, i.e. DLPNO-CCSD(T). [215][216][217] The method has already been applied with success in smaller transition metal systems, 218-221 but the OEC presents much greater challenges. Pilot applications on minimal structural models have been reported but it remains to be seen whether successful applications are possible on realistic models of the OEC.
In the context of electronic structure challenges it should also be noted that the remarks made above regarding the nature of broken symmetry solutions also hold for any other aspect of the system, including DFT-derived mechanistic schemes. The use of a single spin configuration to follow a reaction pathway, particularly a reaction that involves multiple spin states of an exchangecoupled system, is an approximation whose crudeness is   29 Hydrogenases produce H 2 with turnover frequencies reaching (3 Â 10 4 s À1 at 30 1C) at marginal potential beyond equilibrium (À400 mV vs. NHE, pH = 7). 222 This catalytic activity can only be matched by platinum. 223 Their catalytic efficiency, thermostability, and low overpotentials, make hydrogenases ideal models for the bioinspired synthesis of new organometallic catalysts. However, their complex protein structure and the nature of the reaction they perform have posed obstacles in obtaining a clear and precise view on the details of how they function. The size and composition of the active sites make hydrogenases challenging systems to study also from the theoretical point of view. 224 The structural elucidation of the [NiFe] hydrogenase revealed an organometallic Ni-Fe complex with CO and CN-ligands bound to Fe, 225 while four cysteine residues ligate the Ni site. Two of the cysteines bridge the two metal ions as shown in Scheme 1b. The nickel center has been suggested to play a critical role in hydrogen binding/proton activation and in the redox chemistry of the catalytic cycle. It is postulated that three different active states are involved during catalysis: Ni-SI a , Ni-R, and Ni-C. 34 Although broad agreement exists on the presence of Ni-C and Ni-R, uncertainties remain about the completion of the cycle and about how each state converts to the other.
Consequently, there is still lack of consensus on the catalytic cycle of [NiFe]-hydrogenases even if distinct mechanisms have been proposed on the basis of combined experimental and theoretical evidence (Fig. 10). In the proposed cycle, H 2 is activated through binding to the diamagnetic Ni-SI a state, where the metals are present as Ni II and Fe II and have a low spin state. This activation leads to heterolytic cleavage of the H-H bond and the formation of the Ni-R state, which contains a bridging hydrido ligand. The two-electron oxidation of Ni-R regenerates the Ni-SI a while the second proton is released. The intermediate Ni-C species corresponds to the one-electron oxidized form of Ni-R and is assigned to a Ni III Fe II -hydrido complex, which has been experimentally characterized. 226 The  227 The ligand sphere leads to Fe II metal centers in the low-spin state. A key structural difference to the [NiFe]hydrogenase is the presence of a dithiomethylamine (dtma) ligand that bridges the Fe ions. 228 The catalytic cycle of [FeFe] HG (Fig. 11) remains controversial despite extensive research efforts. According to recent reports, binding of H 2 to the oxidized active ready state, H ox , results in the heterolytic cleavage of H 2 , and in the formation of H hyd . The latter species is a mixed-valence Fe II Fe I intermediate featuring a protonated dtma and a terminal hydride ligand. On deprotonation, H hyd is converted into the rather stable H sred state with both metals being at the (+I) oxidation state. The one electron   oxidation of H sred yields H red before a proton-coupled electron transfer recycle H ox . Overall, the mechanism of the [FeFe] HG remains debated particularly with respect to the one-and twoelectron reduced intermediates. 33 The absence of consensus regarding the complete catalytic cycle of hydrogenases has driven developments towards proficient catalysts for HER. Studying smaller catalysts is insightful for the proper design of new molecular electrocatalysts even if they lack the functionally important protein matrix. Many biomimetic and bio-inspired complexes have been developed and investigated with the purpose to understand and mimic the proficiency of hydrogenases. In this respect, computational approaches combined with experiment can play a key role in the development of efficient synthetic catalysts. [229][230][231] Contributions of quantum chemistry are not limited to rationalizing experimental observations; the correlation between theory and experiment drives conceptual and methodological developments toward new experimentally accessible systems. As in the case of the OEC, studies combining experiment and theory contribute to understanding the properties and reactivity of the enzymes, but they have also assisted more directly in the design of more effective HER catalysts.
In the following, we report on the extensive use of the theoretical methods to better apprehend and get insight into the geometric structures, electronic structures, and properties of hydrogenases as well as related biomimetic and bio-inspired models. We will show that quantum chemical studies are crucial to identify the electronic and thermodynamic parameters that govern the reactivity of catalytic centers, to determine the crucial structural determinants of catalytic activity, and to unravel the key steps of their reaction mechanism. Such investigations in the fields of hydrogen evolution nowadays offer ample proof of the use of quantum chemistry as a predictive tool and as an essential component in the understanding of the natural system and in the design of novel synthetic catalysts.

Geometric structures
The complexity of hydrogenases and the various limitations pertaining to their experimental characterization have fuelled great interest in applying theoretical methods to investigate the structures of models of the different redox states of HG. Looking at the structural aspect of the hydrogenases, theoretical tools have been successfully employed to model both the models and the enzymes. It is now accepted that geometries can be reliably predicted by DFT methods and converge quickly with basis set size, making such prediction rather economical. Long experience has established that the cluster approach using relatively small models (for example, up to ca. 200 atoms) of the protein active sites can provide meaningful mechanistic insights like hydrogen evolution by HG (Fig. 12). 53,54 On the other hand, realistic modelling of the protein environment is achieved by quantum mechanics/molecular mechanics (QM/MM) and ''Big QM'' approaches. In the case of the hydrogenases, many studies employing these approaches have been reported over the years and contributed to understanding the enzymatic environment and the reaction channels of hydrogenases. 226,232 In the case of Computational studies have identified the heterolytic cleavage of H 2 as the initial step of the reaction and determined that one hydrogen becomes a bridging hydride whereas the other is accommodated as a proton on one of the cysteine residues. 233,234 For [FeFe] hydrogenase, QM/MM studies have been extensive and multifaceted. For example, it was shown that the inclusion of two [4Fe4S] clusters in the QM part in addition to the H-cluster provided useful insights into the oxidation states of the active site and the interplay between the two cubanes. [235][236][237] We can also cite alternative approaches such as the quantum refinement approach based on accurate DFT calculations that was applied to the [FeFe] enzyme to elucidate the nature of the dithiolate ligand. 238 Molecular dynamics studies have offered support for a proton transport channel to facilitate substrate and product delivery in both hydrogenases. 239 A specific pathway was recently suggested for the second proton, whereby reduction of the active site triggers a structural change, opening a water channel with favourable energetics for the second proton transfer in the [FeFe] HG. 240 In addition, a glutamic residue near the active site was shown to be important for proton transport as its mutation results in loss of activity of [NiFe] hydrogenase. 241 Molecular dynamics can be also combined to other approaches like the coarse-grained (CG) analysis in order to probe important catalytic steps of [FeFe]-hydrogenase, including the transport of electrons to the active site via the iron-sulfur clusters. 242 Population analysis schemes are simple yet valuable tools, as recently proven by a study using conceptual DFT and Natural Bond Order (NBO) to get new insight regarding the transformation of a proton into a hydride during the catalytic cycle of the [NiFe] HG. Their results indicated that such a transformation is driven by spontaneous rearrangements of the electron density at room temperature. 243 Such analysis can be very informative about the ligand binding mode and strength to the metal centers of the [FeFe] active site. For instance, a recent report showed that, upon reduction, the distal iron is favoured by the bridging carbonyl which suggests that electrostatic stabilization would influence the charge accumulation occurring prior to H 2 release. 244 Topological analyses based on NBO were also used to revisit the nature and strength of the hydride and dihydrogen interaction with the iron center in the [Ni-Fe] enzyme. 245

Electronic structures
Deciphering the spin states of the species that are involved in the catalytic cycle of hydrogenases has reached a whole new level with recent studies showing that it is now possible to conduct calculations with accurate wave function quantum chemical methods, such as the coupled cluster method with single, double and perturbative triple excitations, CCSD(T), and multireference approaches (CASSCF, DMRG). These new advances are very important as energies for reactions involving two-electron transfers to the metal centers like in hydrogenases can be sensitive to the nature of the DFT functional and higherlevel methods are thus necessary to obtain more reliable energies and reaction pathways. For example, multi-reference Møller-Plesset 246 and CAS-PT2 calculations 247 were conducted on the active site of the [NiFe] hydrogenase and helped to assign the ground spin state as a singlet through the evaluation of singlet-triplet energy splitting. Wave function methods including CCSD(T) and DMRG have been employed to investigate the interaction of dihydrogen with the [NiFe] active site. These calculations showed that H 2 binds more favourably the Ni site in the singlet state and provided additional understanding about how the active site behaves in the presence of H 2 . 248 A similar work based on multiconfigurational short-range DFT calculations supported this conclusion. 249 DMRG has also been applied to study the electronic structure of [4Fe-4S] clusters relevant to the description of the hydrogenase active sites 125 and demonstrated the importance of developing sophisticated methods for electronic structure calculations for an accurate evaluation of the low-lying states that have implications for the reactivity of most bioinorganic systems.
A different aspect of spin-state energetics is found with magnetically interacting (exchange-coupled) transition metal ions as in the case of the [FeFe] hydrogenase whose active site, the H-cluster, corresponds to a [2Fe] subunit coupled to a [4Fe4S] cluster. Such systems can be handled within the DFT framework using the broken-symmetry (BS) approach, which gives access to exchange coupling constants and hence magnetic spin states and energies. Several studies have been reported using this approach showing its capability to provide useful insight on the electronic configuration of the complicated [FeFe] active site by revealing the distribution of metal valences and spin-coupling schemes for the iron-sulfur clusters. [250][251][252] Similar studies were conducted in the case of the O 2 -tolerant [NiFe] hydrogenases that contain a proximal [4Fe-3S] cluster whose redox states and magnetic couplings were elucidated and discussed using BS-DFT calculations. [253][254][255]

Spectroscopic properties
In conjunction with experimental methods, theoretical calculations appear as an invaluable tool for the interpretation of spectroscopic results. Spectroscopic data for a sufficiently complex system rarely have a unique interpretation in terms of electronic structure. On the other hand, it can be difficult or impossible to identify the most appropriate electronic structure description using only total energies as criterion in quantum chemical studies. Therefore, the explicit calculation of spectroscopic properties is the way toward a better understanding and more reliable verification of the geometric and electronic structures of a complex system that is being modeled. Spectroscopic techniques including FT-IR, Nuclear resonance vibrational spectroscopy (NRVS), Electron Paramagnetic Resonance (EPR) or Mössbauer have been widely employed to investigate [NiFe] and [FeFe] enzymes and models, and often combined with theoretical methods (Fig. 13). 29 Following the development of a computational methodology for an accurate prediction of the IR spectra of the active site of [FeFe]-hydrogenase, 256 many studies have been conducted using DFT calculations and helped in the unambiguous assignment of the structure and redox states of intermediates like for the H-cluster. [257][258][259][260] While the calculated absolute vibrational frequencies are often not quantitatively accurate, shifts in frequencies due to protonation and/or reduction events are reliable. Due to the peculiar vibrational frequencies of the CN and CO groups, it is possible to distinguish the redox states of the [NiFe] active site. Indeed, the stretching frequencies of these ligands are spectroscopic fingerprints that can probe the electronic structure of the metal core and the structural composition of the active site. 261 Comparison between computed and experimental vibrational frequencies can provide even more subtle information like the influence of the second coordination sphere on the electronic structure of the [NiFe] active site. 262 A review article recently highlighted that the correlation between experimental and DFT-calculated IR spectra is a powerful tool to identify some key features of the active site cofactor of HG and to unravel the structure of transient intermediates with great precision. 227 NRVS is emerging as spectroscopic technique to study hydrogenases. The vibrational density of states spectra that are obtained when using NRVS can be reliably predicted using density functional theory (DFT) and assists in spectral assignments. A striking example comes from a recent study that combined IR and NRVS spectroscopic techniques to DFT calculations to elucidate the structures and the local coordination spheres of the one and two-electron reduced intermediates of the [FeFe] HG. Their results showed no evidence for a bridging hydride which may have important consequences for the mechanism of H 2 conversion by this enzyme. 263 266,267 while a recent report showed the direct evidence for a hydride bridge in the active site of the fully reduced Ni-R form whose electronic structure could be assigned as a low-spin Ni 8 (m-H)Fe 8 core. 268 EPR spectroscopic techniques have been demonstrated as extremely useful to investigate paramagnetic intermediate states during hydrogen conversion catalytic cycle. 228 Understanding the hydrogenase mechanism implies being able to determine the exact structure of the active site and its surrounding for all states. In addition to this information, EPR spectroscopy can also give access to the oxidation and spin states of the metal ions involved and their interaction with each other or with the ligand sphere. One known example of the use of theoretical EPR spectroscopy to study HG deals with the Ni-L state, detected when the Ni-C state is put under illumination. DFT calculations were conducted on all possible structures of Ni-L state and g tensors together with hyperfine coupling constants were computed as indicators to compare with EPR data. From the calculated spectroscopic parameters, it was proposed that the proton most likely binds to the terminal cysteines rather than the bridging ones and that a metalmetal bond between nickel and iron was present. 269 Another example lies in the study of an artificial [FeFe] hydrogenase for which DFT calculations helped to determine the structure of the active site. Specially, computation of the hyperfine coupling constants reproduced experimental data and highlighted the presence of a CN ligand with the cyanide C atom bridging one mixed-valent iron of the [4Fe4S] cluster and the N atom bound to one iron of the binuclear unit. 270 Several experimental and theoretical investigations of the H-cluster using advanced pulsed EPR techniques and DFT calculations have been carried out. They showed that spin polarization in the dtma ligand can be linked to the asymmetric coordination of the distal iron site with its terminal CN À and CO ligands which would have implication on stabilizing the iron-hydride intermediate in the catalytic cycle. 271 More recently, EPR data together with DFT calculations clearly showed that a specific conformation of the external CO ligand is favoured so that it binds in an apical position on the iron center. These combined studies highlighted the rigid configuration of the CO and CN ligands which is stabilized by H-bonds and has implication on the reactivity of the enzyme as it ensures minimal reorganization energy during catalysis. 272 As an indispensable analytical tool in iron coordination chemistry, Mössbauer spectroscopy is a valuable spectroscopic technique to study hydrogenases. It can probe the spin state of the metal centers in the active sites as well as their local environment like the coordination geometry. For instance, theoretical Mössbauer spectroscopy has been used to investigate the [4Fe3S] cluster proximal to the active site of the O 2 tolerant [NiFe] HG. The calculated isomer shifts and quadrupole splitting helped in discriminating between different interpretation of the data and a clear assignment of both the spin coupling scheme and the redox states within the iron-sulphur cluster was made possible. 273 Combined with theoretical calculations, this spectroscopic technique has been also employed to study diverse [FeFe] hydrogenase intermediates. The putative structures of these intermediates were obtained from QM/MM calculations and the comparison between calculated and experimental isomer shifts could be used to discriminate between different intermediates following a proper calibration based on iron-based model systems. 274 The same approach was also used for studying the coordination environment of an active site model complex. Based on the calculated isomer shifts and quadrupole splitting, they could determine the protonation states, the presence of a bound water and the redox state of the iron center providing meaningful insight regarding the proton transfer mechanism. 275 A more recent study made use of several techniques including Mössbauer spectroscopy and DFT calculations to set-up structure-property relationships and to show how that the reactivity of all HG enzymes could be mimicked using isomers of a [NiFe] bio-inspired model. 276 With the recent development of dedicated quantum chemistry program packages 277-280 the use of theoretical spectroscopy has been democratized with the aim to connect the results of computations to experimental observables since spectroscopic properties remain a precise and sensitive fingerprint of the electronic structure for a given system. Consequently, combining spectroscopic characterizations with theoretical calculations has been successfully employed by many research groups to identify the geometric and electronic structures of biomimetic 40,60,265,268,[281][282][283] and bio-inspired 284-290 models of hydrogenases and appears as a prerequisite to any further computational investigation of their reactivity towards H 2 evolution.

Thermodynamics
Studying the reactivity of hydrogenases and related model systems implies the initial determination of thermodynamic properties like redox potentials and pK a values. Indeed, the hydrogen evolution reaction proceeds via the transfer of two electrons and two protons to form H 2 . For a given catalyst, six sequences can be formally considered: ECEC, ECCE, EECC, CCEE, CEEC and CECE with E being the electrochemical event i.e. the electron transfer step and C being the chemical event i.e. the proton transfer step (Fig. 14). Through the computational evaluation of the free energies of the species associated to these chemical and electrochemical steps, it is possible to get access to the thermodynamic quantities with a reasonably high level of accuracy. The theoretical prediction of these parameters is usually conducted using Born-Haber cycle and the calculated Gibbs free energy change associated with the reduction and protonation event in solution, respectively. To calculate experimentally relevant electrochemical quantities, reference systems have to be considered. To do so, the absolute reduction potential of the reference electrode and the pK a of the acid serving as the proton source have to be computed and further subtracted from the Gibbs free energy change in solution for reduction and protonation. 229,291 Electronic structure methods have extensively used to predict relative redox potentials and pK a values of molecular electrocatalysts and hydrogenases models as a first step for further mechanistic investigations. 33,226,231,287,292,293 One useful tool for envisioning the relationship between E and C steps during hydrogen evolution finally lies in the construction of Pourbaix diagrams. These diagrams depict the thermodynamically most stable species for a given reduction potential and pH value (Fig. 15). For a given catalyst, constructing a theoretical Pourbaix diagram allows to evaluate the stability of all possible intermediate species on a broad pH range. Linear correlations using Pourbaix diagrams, redox potentials and pK a values have been established for cobalt and nickel electrocatalysts as powerful tools in catalyst design, enabling the computational prediction of properties for electrocatalysts that have not yet been synthesized. 292,294-296

Reaction mechanisms
Computational approaches have been widely used to elucidate various mechanistic aspects of catalytic hydrogen evolution, including possible protonation sites, energetics of distinct pathways, order of steps, and stepwise versus concerted   processes. The reliable calculation of relative stabilities of various isomeric forms and the detailed insight into plausible pathways for their interconversion are crucial to address reactivity problems in hydrogenase chemistry. Using DFT methods, free energy profiles of possible mechanistic routes can be evaluated along with the corresponding intermediates and transition states resulting from reduction and/or protonation events. The theoretical determination of the quantities associated to the electrochemical and chemical steps was shown to help to determine which pathway is favoured among all possible sequences. Many studies made use of detailed DFT calculations to inspect free energy landscapes and to identify the predominant pathway for proton reduction among all possible sequences in the case of molecular electrocatalysts 61,62,292,294,[297][298][299][300][301][302][303][304][305][306][307][308][309] and hydrogenase models (Fig. 16). 39,59,281,282,287,288,[310][311][312][313][314][315][316] Environmental effects play an important role in modulating the energetics of every intermediate state in the reaction pathway of HER catalysts. Specially, the presence of solvent molecules or anions from the acid serving as the proton source can play a crucial role in modulating the kinetic and thermodynamic aspects of the reaction. 317 For instance, the proton source can affect the thermodynamics of catalytic event and a theoretical description of these effects on the pK a value, metal redox states or protonation states are necessary to inspect the feasibility of different pathways for the catalytic cycle of cobalt catalysts. 293,318 The DFT investigation of a nickel catalyst reported the effect of the solvent on redox properties and how it can affect the stability of the metal oxidation state by changing its coordination. In this study, accurate estimates of redox potentials could be estimated by accounting for a possible ligand release. 296 The importance of the anation was demonstrated in a recent DFT study on a cobalt catalyst which showed that substitution of a bound acetonitrile with acetate alters the proton reduction energetics (Fig. 17). 317 Thus, accounting for these effects and the possible direct involvement of the solvent are critical considerations in computational studies.
For a complete rationalization of the reactivity towards hydrogen evolution, it is often sufficient to evaluate the activation barrier for thermally activated processes. For other cases, free energies of reactions have to be combined with kinetic rates (Fig. 18). To do so, the rate constants for proton transfer, electron transfer, and concerted events have to be calculated as the kinetics of these various steps can play an important role in the reaction pathways. Such an estimate requires further analysis, involving accurate evaluation of reorganization energies, adiabatic couplings for electron transfer, and methods that can properly sample the quantum dynamics of proton-coupled electron transfer. Several studies on molecular electrocatalysts have reported the theoretical evaluation of electrochemical non-    adiabatic rate constants for electron and proton transfer within the Marcus theory framework. 229,[319][320][321] In the case of hydrogenases, kinetic rates have proven hard to extract. 322 However, combining MD data with phenomenological rate equations, hydrogen diffusion constants in HG were recently calculated and found to be in agreement with experiments. 323,324 Computational studies on the mechanistic aspects of hydrogen evolution by hydrogenases have been an active topic of research in many groups. 224,291,325,326 For instance, various computational works have been performed to study the reaction cycle of the [NiFe] HG with a special emphasis on H 2 binding site and proton binding position. 233,234,327,328 While there has not being a complete consensus on the catalytic cycle (Fig. 10), the mechanism presumably starts with the Ni-SI a state binding H 2 on the Ni II site. H 2 is heterolytically cleaved yielding the Ni-R state displaying a hydride bound to both two metals and a proton binding the Cys546. The next step involves two possibilities, depending on whether a proton or an electron is removed first. Both scenario lead to the formation of the Ni-C intermediate featuring a Ni III center. Subsequently, the formed hydride is oxidized, the metal reaches its Ni I state, and the resulting proton is transferred to Cys546. The cycle is then completed with the formation of the Ni-SI a resting state. The reaction mechanism of H 2 oxidation by [NiFe] hydrogenase was recently revised using DFT methods and focused on the energetics of the entire catalytic cycle. Here two different mechanisms were investigated: the heterolytic cleavage of H 2 and the homolytic one that occur with different redox states of the [NiFe] core. Their main finding is that reaching the homolytic mechanism implies one cycle of the heterolytic mechanism. They also correlated the design of the [NiFe] active site complex to the criteria of energy minimization upon proton and electron release, and thus observed an endergonic binding of H 2 which has direct implications on the development of efficient HG mimics. 329,330 Another study using DFT calculations recently reported new mechanistic insights on the reactivation process involving oxidized forms of the [NiFe] HG. While [FeFe]-hydrogenases are irreversibly inactivated by O 2 , [NiFe] are more resistant to oxidation. However, oxygen exposure creates a mixture of two inactive, oxidized forms, the Ni-B and Ni-A states. Importantly, these are not permanently inactivated but can be reactivated by reduction. The results of this study showed that the bridging hydroxide ligand found in the two states is lost upon reduction and protonation steps at the active site. Structures for the species obtained upon reduction of Ni-A and Ni-B were proposed and specific structural features were shown to impact the reactivation kinetics of the oxidized and inactive states. 331 Similarly, with [NiFe] hydrogenases, distinct mechanisms were proposed for the [FeFe] ones and differ mainly in the protonation sites and H 2 binding. Lately, computational studies together with experimental investigations seem consistent with one scenario. When considering the catalytic cycle in the direction of H 2 production (Fig. 11), the [FeFe]-hydrogenase state directly involved in hydrogen binding is H ox . In this state, the redox of the H-cluster can formally be described as Fe II Fe I . One electron reduction and protonation of H ox yields a transient Fe I Fe I species featuring a protonated dtma ligand and referred as H red . Further proton-electron transfer takes place resulting in the formation of H hyd , a mixed-valence Fe I Fe II species featuring a terminal hydride ligand and a protonated amine. H 2 detachment from this intermediate brings the H-cluster back to the initial Fe II Fe I redox state. [332][333][334][335] The catalytic activation of H 2 by [FeFe] HG has been recently revised by means of DFT methods to elucidate the specific details underlying the different performances of the enzymes when compared to biomimetic models. Their results highlighted that the enzyme performs better which can be correlated with two requisites: (i) the presence of an electronpoor iron center which favors H 2 binding and (ii) a electron-rich cofactor which ensures oxidation of the metal-H 2 adduct. Disclosure of these differences will guide the development of novel biomimetic models and lead to reconsider the utility of CN ligand in the catalyst design. 336 In addition to the extensive use of DFT methods, QM/MM thermodynamic cycle perturbation (QTCP) computations were recently shown as a valuable approach to obtain accurate free energies to model the catalytic cycle for hydrogen evolution by hydrogenases. In these calculations, the QM/MM method is used to compute the energy differences on configurations generated by MD simulations. QTCP currently gives access to reaction energies with an accuracy of 20 kJ mol À1 . Reaching higher degree of precision will require, among other things, to take into account other factors like polarization or short-range electrostatics and the development of high-level methods, both being the topic of intense research efforts. 233,234,337,338

Conclusions
Impressive advances have been made in recent years in understanding the structure and function of the photosynthetic oxygen-evolving complex and of hydrogenases. Quantum chemistry has played a key role in most of these efforts, especially when it has been fruitfully combined with experiment. Nevertheless, photosystem II and the hydrogenases continue to pose great challenges, both in terms of fundamental understanding and in terms of the technical difficulties involved in studying them. In the field of water oxidation an incredible amount of insight has been gained on the geometric and electronic structure of the oxygen-evolving complex, yet serious questions remain open regarding the precise nature of several catalytic intermediates and nothing is known with certainty about the critical dioxygen formation and evolution steps. In the field of hydrogen evolution, despite the recent advances made in the understanding of the mechanism by which hydrogenases exert their function, especially the formation of hydride species or the key role of the metal site involved in catalysis, many essential aspects of the catalytic process including CO/O 2 inhibition or activation/deactivation processes are incompletely understood. From the point of view of quantum chemistry both more accurate methods and more reliable models will be needed to retain the pace of theory-led discovery in the field of metalloenzyme solar fuel research. Among the points highlighted in this article is the development of first principles electronic structure methods that overcome deficiencies associated with existing DFT implementations and deliver accurate predictions of spin-state and reaction energetics for systems containing several interacting open-shell metal ions. Coupled to that will be the more extensive and more reliable use of quantum chemical approaches for the prediction of spectroscopic properties, which is a major point of contact between electronic structure theory and experiment. In terms of model definition, the conventional QM-cluster approach will be increasingly disfavoured over a more complete treatment of the protein environment. This is not a luxury but a necessity: the active site is not a catalyst, the whole enzyme is. The protein matrix is exclusively responsible for critical elements of enzyme function such as electron transfer, regulation of substrate delivery, proton access/egress, product release, etc., therefore our understanding of enzyme function will keep expanding to encompass the functional role of the protein beyond the immediate vicinity of the active sites. A multitude of different methods are emerging to achieve this type of treatment, many of which will in addition be able to offer access to a dynamical view of the system. Computational simulations will increasingly influence our understanding of kinetics related both to the transport of substrates to active sites that are buried within the protein and only accessible via specific channels, and to the release of products. With the broader goal in mind of overcoming the efficiency gap between the biological systems and biomimetic/bio-inspired compounds, the most successful strategy remains to combine insights from both systems to uncover differences as well as common principles that can be put to use in the design of water oxidation and hydrogen evolution catalysts. Computational methods have made tremendous advances to reach the level of utility they have today in this effort, and their contribution in the future is certain to increase.

Conflicts of interest
There are no conflicts to declare.