How ab initio Molecular Dynamics Can Change the Understanding on Transition Metal Catalysed Water Oxidation.

Artificial water splitting is a promising technology that allows the storage of renewable energy in the form of energy-rich compounds. This mini-review showcases how theoretical studies contribute to the under-standing of existing water oxidation catalysts (WOCs) as well as inspiring the development of novel WOCs. In order to understand the chemical complexity of transition metal complexes and their interaction with the solvent environment, the use of sophisticated simulation protocols is necessary. As an illustration, a family of ruthe- nium-based WOCs is presented which were investigated employing a wide range of forefront computational methods with emphasis on ab initiomolecular dynamic based approaches. In those studies a base assisted oxygen-oxygen bond formation was identified as the energetically most favourable reaction mechanism. By examining the role of local environmental effects at ambient temperature and the effect of modifications in the ligand framework, a comprehensible picture of the WOCs can be given, where the latter can serve as a guideline for further experimental and computational studies. In this mini-review, we provide a description of the methods, and the findings of our previous computational studies in compacted form, aimed at scientists with a theoretical as well as experimental background.


Introduction
Artificial water splitting is among the most promising approaches to store renewable energy in the form of chemical bonds. The energy storage in the form of high energy chemical compounds has the possible advantage that those products could be transported conveniently over large distances. [1,2] This would resolve some of the problems of many renewable energy sources, i.e. their geographical and seasonal ties. The water splitting reaction can be divided into two half-reactions: On the one hand, there is the hydrogen production, where protons are reduced to molecular hydrogen. [3] This is usually achieved via the formation of metal hydride species, that upon protonation release the desired product. On the other hand, there is the water oxidation which involves stripping oxygen off its additional electrons which are eventually used in the hydrogen reduction reaction. Due to the harsh condi-chemical transformation. In the past years considerable progress has been made to determine the relevant geometrical descriptors by using machine learning, however those approaches are beyond the scope of this mini-review. [11] CVs are not limited to geometrical parameters, but other chemical descriptors, such as the occupation of the frontier orbitals, can be used to describe the progress of the reaction. [12] Choosing an appropriate low dimensional CV-space is the most crucial step when using most of the existing enhanced sampling methods. The goal is to reduce the computational cost by only sampling a selected region of the phase space. In the following, two of these enhanced sampling methods are briefly introduced.

Bluemoon Ensemble
In 1989, Carter et al. introduced a method, which is referred to as Bluemoon (BM) sampling, which allows simulating rare events in the condensed phase. [13] The theoretical background of BM relies on the scheme of thermodynamic integration to reconstruct the free energy profile of the chemical transformation (∆F). This is achieved by first discretising the CV ([ξ 0 ...ξ 1 ]) for the transformation from reactant to product and then running independent constrained MD simulations for each intermediary value of ξ in order to obtain the Lagrange multiplier (λ). The latter corresponds to the force acting on the system satisfying the given constraint. The numerical evaluation of this free energy difference can be expressed as: where T is the temperature, k B is the Boltzmann constant, and λ defines the strength of the constraint force propagated at each simulation step. The remaining two terms Z and G are given by: where m i refers to the atomic mass of the i th nuclei and r i is the corresponding vector of atomic position. For example, if the reaction is a simple cleavage or formation of a chemical bond, then it is common to use the interatomic distance between atoms i and j as the CV. For this special case, Z and G are independent of r i . [14]

Metadynamics
In this section we describe an enhanced sampling technique called metadynamics (metaD) which was introduced by Laio and Parrinello in 2002. [15] The underlying principle of metaD is to sample the unexplored areas of the phase space by depositing a history dependent bias potential in areas that were already sampled. For reasonably long simulation times (ξ, t → ∞), the history dependent bias potential V(ξ, t) is guaranteed to converge to the free energy of the system with an opposite sign. [15] Conventional metaD suffers from convergence issues related to overfilling of local minima which leads to a highly corrugated FES. The well-tempered metadynamics (WT-metaD) formalism overcomes this issue by adjusting the bias potential deposition stride over the course of the reaction, ranging from a high deposition frequency at the beginning of the simulation to a lower frequency for a nearly converged system. The same effect can also be achieved by reducing the height of the bias potential over the course of the simulation. [17] (1) tions required for this step, this half-reaction is considered as the bottleneck of the whole water splitting process (see Scheme 1) Although there are numerous heterogeneous and homogenous catalysts available, none of them is reported to have the efficiency, long-term stability, and cost-efficiency needed to be economically feasible on a global scale. [2,4] In order to reach this goal it is necessary to understand the thermodynamics, kinetics, and mechanisms governing those reactions at an atomistic level. The spatial and temporal limitations of existing experimental techniques make the characterization of illusive catalytic intermediates a challenge. Computational modelling is thus a valuable tool to bridge this gap and provide a molecular perspective of the water splitting reaction. Nowadays, electronic structure calculations based on wavefunctions and electronic densities are routinely used to supplement the experimental findings. However, these methods are usually used in a static approach, thus not taking into account finite temperature effects, dynamic solute-solvent interactions and ensemble averaging. [5][6][7] When it comes to ambient temperature effects and solute-solvent interactions, first-principles dynamic methods, which go beyond those standard static approaches, are desirable for a more realistic modelling of the system of interest.
In this short review we will present some recent computational findings in the field of homogeneous transition metal-based water oxidation with emphasis on recent work in our group relying on density functional theory (DFT)-based molecular dynamics (MD), often termed ab initio MD. In the following section we will briefly introduce some of the common simulation protocols for DFT-MD which is followed by the discussion of the application of these methods to the water oxidation reaction.

Methods
Chemical transformations where an equilibrated system transits from the reactant to the product state, passing through a high energy transition state, are considered to be rare events in the context of MD simulations. Therefore, in order to accurately describe the underlying free energy surfaces (FESs) suitable simulation protocols are required. [8] Regarding the electronic structure methods, a number of wavefunction-based methods have become invaluable for working toward understanding reaction mechanisms at the atomistic scale. [9] However, the complexity of these methods and involved computational cost, particularly in the context of large systems, have limited their application to MDs. On the other hand, DFT offers to be a reasonable compromise between computational efficiency and accuracy. [10]

Collective Variables
Chemical transformations evolve usually only in a small subspace of the 3N-dimensional phase space (N corresponds to the number of atoms in the system). Computationally, the term used for a variable describing a reaction coordinate capturing the progress of the reaction in a subspace of the whole 3N-dimensional space is collective variable (CV). Most often, CVs are intuitively chosen as geometrical descriptors such as distances or angles between reacting molecules and/or atoms, that best describe the Scheme 1. Water splitting half-reactions and their standard one-electron reduction potential (E°) (vs. NHE, pH=0).
ifications in the ligand environment were proposed that possibly lead to a better catalytic performance. We have used a similar approach previously to study a model Co-based WOC. [30]

Ligand Acidity
Acid dissociation constant (pK A ) and reduction potential are important parameters which can be used to elucidate the behaviour of catalysts. Nevertheless, for most intermediates of the catalytic cycle they are experimentally inaccessible. Despite the fact that there are numerous simulation protocols available, [31,32] calculation of those properties still remains a challenge. Particularly, since in many cases the use of an implicit solvent model is insufficient to achieve chemical accuracy. An atomistic and dynamical description of the solvent in a simulation box obeying periodic boundary conditions is therefore attractive although it comes with a significantly higher computational cost due to the increase in system size and the requirement for statistical sampling of the configurational space.
Some of the most popular enhanced sampling protocols for calculating dissociation constants are based on a thermodynamic integration scheme. On one hand, there is the 'insertion-deletion' (I-D) technique which was developed by Sprik and co-workers. [33] The main idea of this method is to calculate the free energy associated with the deprotonation of the acidic site (∆ pt F AH ) by removing the proton from the simulation cell (deletion step). Similarly, in the insertion step, the free energy of the protonation of a solvent molecule (∆ H20+ ), i.e. the formation of a hydronium ∆F H30+ , is calculated. The pK A values are then obtained by combining those two free energies according to Eqn. (6): The I-D technique was successfully applied to various systems, [34] most recently by Govindarajan et al. to Ru-based WOCs. [35] Analogous simulation protocols can be used to calculate reduction potentials or dehydrogenation free energies, which has been carried out for a cobalt aquo ion in the context of water oxidation. [36] On the other hand, there is the constraint DFT-MD (BM) based approach. The basic principle of this method is described in section 2. Here we only discuss its application to determine pK A values. In literature there are alternative views on how exactly the free energy difference (∆F) obtained by this method is related to the pK A . The proposed post-processing protocols mainly differ from each other in terms of the definition of the reference reaction, i.e. the protonation of a different basic or solvent molecule (the reader is referred to ref. [28] for more details). However, their performance has never been directly compared, in particular not in the context of transition metal complexes.
We have determined pK A values not only for small molecules such as formic acid and phenol but also for derivatives of RuPy5 WOCs by using different post-processing protocols. While all methods were able to qualitatively reproduce experimentally determined pK A values, the probabilistic approach introduced by Davies et al. [32] was found to be the most reliable one. It is worth mentioning that the performance of the pK A BM method deteriorated when it came to predict pK A values of highly acidic species. A possible reason for this lies in the definition of the applied constraint which in our study was the distance between the acidic site and the proton. This choice might not be justified for all cases since proton transfer might happen by the Grotthuss mechanism. [37] To circumvent this problem, different CVs such as the coordination number of the acidic site can be employed. [38] The DFT-MD based pK A studies, where a transition metal complex is modelled with explicit solvent molecules in a periodic simulation box, showcase the fact that nowadays both the simulation protocols and the avail- Within the WT-metaD formalism the expression of the history dependent bias potential is given as The simulation time ξ i (t) is the value of the i th out of d CVs at time t, τ is the potential deposition stride. Furthermore, the Gaussian potential is composed of the following parameters: the initial height (w 0 ) and the width of the Gaussian (σ i ), and an initial temperature ∆T that controls the scaling of the Gaussian height.
BM and metaD have either been implemented directly into various MD program packages or can be used by means of the PLUMED software which has been interfaced with many codes. [18] 3. Ruthenium-based Polypyridyl Complexes

The System
Since the discovery of the 'blue-dimer' (cis,cis-[(bpy) 2 (H 2 O) Ru III ORu III (OH 2 )(bpy) 2 ] 4+ , where byp = 2,2'-bipyridine) in 1982, notably the first molecular artificial water oxidation catalyst (WOC), numerous transition metal-based WOCs have been proposed. [19] Despite the great efforts made to develop WOCs containing cheaper alternatives to noble metals, [20,21] today's most potent WOCs are still based on noble metals such as ruthenium or iridium. [4,21,22] In this context it is worth to mention a family of 2,2'-bipyridine-6,6'-dicarboxylate (bda) based mononuclear Ru-WOCs introduced by Sun and co-workers. [23,24] Derivatives of the latter have been reported to be the most efficient Ru-based WOCs up to date. [24] Besides bda ligands there is a growing number of polypyridyl ligands that were found to stabilise the high oxidation states of the Ru metal centre during the course of water oxidation reaction. [25] Among them is the family of pentapyridyl ligands (py5) introduced by Gil-Sepulcre et al. [26] The ability of Ru-complexes based on these py5 ligands (RuPy5) to oxidize water was later studied thoroughly. [16] Our DFT-based simulations have suggested a mechanistic proposal that is in agreement with experimental findings, see Fig. 1 for an illustration of the proposed active form of the catalyst. [26,27] Further we have shown, by using a rational in silico design approach, the complex interplay between electronic and steric factors that govern the reaction. [27] Based on these findings, mod- Fig. 1. Illustration of RuPy5 and its derivatives. Left: Hydrolysed catalyst that is inactive as a WOC. The pK A of the aquo ligand was calculated with the Bluemoon simulation protocol. [28] Right: Active form of the catalyst. Only two derivatives of the catalyst were synthesized (L 1 = OMe and Me, L 2 , L 3 , L 4 = H). The other derivatives were proposed based on in silico design. [27] The most promising candidate (L 1 = OMe, L 2 = OMe, L 3 , L 4 = H) was then studied further. [29] bilization would require several solvent molecules. Accordingly we investigated the system by using DFT-MD in a simulation box consisting of the catalyst and 107 water molecules (see Fig.  3). [42] Chemical reactions such as the O-O bond formation are rare events within the achievable simulation times and thus, similar to the determination of the pK A values, enhanced sampling methods were employed.
The BM methodology promises to be a reasonable and straight-forward first choice to model the O-O bond formation reaction since chemistry dictates the otherwise difficult choice of a reasonable CV. Employing the O-O distance as the CV we found that the base-assisted WNA reaction mechanism is indeed energetically favoured over the base-independent case (see Fig. 4).
This differences in the mechanism cannot be identified using static calculations to optimize TS structures, due to the formation of a hydronium. In literature there are several examples of BM method applications for modelling the WNA reaction. [40,[43][44][45] A common inference is that the presence of an intra-or intermolec- able computing resources are sufficient to study large and complex homogeneous systems. Further, chemical properties of illusive intermediates can be predicted in order to not only better understand the underlying reaction mechanisms but also to fine tune them to further increase the catalytic performance. This also applies to the oxygen-oxygen (O-O) bond formation, which is a key step in water oxidation and discussed in the next section.

Oxygen-Oxygen Bond Formation
The O-O bond formation is often considered to be the bottleneck in water oxidation catalysis. In-depth understanding of the latter is therefore of great importance to further improve or develop novel WOCs. From a mechanistic point of view, the O-O bond formation is assumed to follow one of two possible pathways. On the one hand, there is the radical coupling (RC) mechanism which involves two metal-oxo species (see Fig. 2, top right corner in pale blue). On the other hand, there is the so-called water nucleophilic attack (WNA) mechanism. [5] Here the metal-oxo species is assumed to undergo a nucleophilic attack by a solvent water molecule(see Fig. 2, left half in blue). The classification of the O-O bond formation into those two mechanisms is rather coarse. At an atomistic level there are many variations of those two basic mechanisms depending on the nature of the ligand framework and the composition of the reaction mixture. For example, the presence of an intra-or an intermolecular base like pyridine can enhance the nucleophilicity of the water molecule by deprotonating it. Pushkar et al. envisioned a different role for dangling pyridyl ligands. It was proposed that these pyridyl moieties can be oxidized by the metal-oxo in an oxygen-atom transfer (OAT) reaction forming a pyridyl-N-oxide (see Fig. 2, bottom right corner in red) that might be involved in the O-O bond formation. [39,40] The variability in terms of possible reaction mechanisms of the O-O bond formation is a challenge for computational chemists attempting to model the reaction, in particular the WNA which involves not only the O-O bond formation but also a proton transfer from the nucleophile to either a basic site or to the solvent. As a consequence, to model the transition state (TS) of the O-O bond formation the inclusion of at least a few explicit solvent molecules is necessary. In conventional TS optimizations, i.e. based on geometry optimizations, this is a valid approach as long as the number of solvent molecules is small. [41] In the first stage of our studies, we used this static approach to model the TSs of RuPy5 WOCs. [16,27] The calculations suggested that the unique ligand framework featuring a dangling pyridyl subunit can act as an intramolecular base, and a base-assisted WNA is feasible. Notwithstanding, the static calculations cannot be used to distinguish a base-assisted and base-independent WNA mechanism, as the latter involves the formation of a hydronium ion, whose sta- ular base facilitates the reaction. [42,43] This can be explained by the fact that the formation of a hydronium is excluded a priori. In the base-assisted WNA not only the O-O bond, but also an N-H bond are formed in case of the investigated RuPy5 WOCs. This implies that a simple CV such as the O-O distance is not reliable to model the reaction in this case. [42] In principle there are no limitations towards the types of CVs that can be employed in BM except for the fact that the Z and G terms (see Eqns (2) and (3)) have to be calculated explicitly for every configuration of the MD trajectory. As a CV we used the difference of the O-O and O-H bond lengths to model the WNA including the proton transfer. However, the fact that the two protons of the nucleophile are indistinguishable at larger distances from the intramolecular base demonstrated a limitation of this simple CV. The BM method can only be used to simulate FESs that depend on a single CV. This limits its application to complex chemical processes, and it is therefore recommended to use methods that can inherently sample the phase space along multiple CVs.
One of these methods is metaD, the methodological details of which are introduced in section 2. While popular in the field of computational chemistry, especially in the context of classical force-field based MD simulations, the method has not seen lots of attention in the context of transition metal-based water oxidation catalysis. Nevertheless, there are a few studies where DFT-MD based metaD was used to study the O-O bond formation facilitated by mononuclear as well as tetranuclear Ru-based catalysts. [44,46,47] In these studies, metaD was used to find a low energy transition state connecting the reactant to the product state. This transition path was further refined using methods such as BM [44] or single point electronic energy structure calculations using a hybrid exchange-correlation density functional. [47] In our recent work we have gone beyond those approaches and modelled the O-O bond formation by a base-assisted WNA with the aim of exploring the minimum energy path (MEP) and the entire phase space relevant for describing the studied reaction. We have used two CVs to monitor the O-O bond formation, the protonation state of the nucleophile and of the intramolecular base. The latter was achieved by using the difference of two coordination numbers as our second CV. The first coordination number describes the coordination of the relevant protons to the nucleophile, i.e. the number of protons bound to the nucleophile, while the second coordination number monitors the protonation state of the intramolecular base. This combination of CVs allowed us to follow the reaction progress and solve the ambiguity concerning the two protons discussed before.
The FES reconstructed from the metaD simulation using those two CVs was further analysed in detail. This has led to a detailed description of the base-assisted WNA and the identification of the MEP, and gave valuable information on the reactivity of the product state. The hydroperoxo moiety formed during the WNA was found to be prone to undergo proton-transfer reactions assisted by the intramolecular base and the solvent. This is in accordance with the commonly assumed next steps of the water oxidation reaction (see Fig. 2), which are a deprotonation reaction followed by proton-coupled-electron-transfer prior to the release of molecular oxygen.
In one of our previous studies, we proposed some modifications in the Py5 ligand framework which according to static simulation has a lower activation barrier for the WNA as compared to the WOC. [27] The most promising candidate was found to have a methoxy group in the para position of the dangling pyridine (Py5OMe (see Fig. 1, Py5: L 1 = OMe, L 2 = OMe, L 3 , L 4 =H). This modification supposedly increases the basicity of the intramolecular base which in turn could facilitate the deprotonation of the nucleophile. [48] Employing the same metaD based simulation protocol as previously in the case of RuPy5, [42] we rationalized how the increased basicity of the dangling pyridine affects the MEP. [29] At first glance, the FESs of RuPy5 and RuPy5OMe were found to be reasonably similar, i.e. the same local minima were explored, and so were the connecting MEPs (see Fig. 5). However, by projecting the MEP onto a single CV, i.e. the O-O distance, it became possible to identify three distinctive stages of the reaction (see ref. [29]). In the beginning there is the pre-organization where the nucleophile is sterically orientated for the O-O bond formation and the proton transfer reaction. Secondly, there is a partial proton transfer, where the proton is shared between the nucleophile and the intramolecular base, and finally there is the O-O bond formation. The biggest difference in terms of free energy for the two catalysts was found for the first state of the reaction. In the case of RuPy5OMe this energy was found to be about half of the one obtained for RuPy5 (15±1~kJ/mol vs. 34±1~kJ/mol). [29] Careful evaluation of the FES revealed that this decrease in free energy is caused by the fact that the increased basicity favours solvent configurations where the nucleophile remains in proximity to the intramolecular base by forming a hydrogen bond. In the case of RuPy5OMe we found the associated free energy differences for the later two phases of the reaction also to be lower (3-6~kJ/mol). The fact that the solvent plays an important role in stabilizing the TS explains why in the case of the static TS calculations only marginal structural differences among the two WOCs were found. [27] A long-standing question related to the O-O bond formation is whether the reaction follows a simultaneous two-electron step or consists of two sequential one-electron transfer steps. [49] The important difference between the two mechanisms is the (potentially short-lived) existence of a H 2 O •+ or a hydroxyl (OH • ) in proximity to the metal oxo species. These radicals are supposedly formed by an electron transfer from the nucleophile to the metal-oxo species prior to the actual O-O bond formation. Because single determinate DFT lacks the ability to accurately describe multideterminational electronic states, [7] an approach to analyse entire reaction pathways by means of complete-active-space-self-consistent-field (CASSCF) was proposed recently and applied to follow the electronic configurations during the O-O bond formation of the RuPy5 WOC. [50] These calculations suggested that both one-and two-electron processes are possible. Nevertheless, the mechanism is dominated by a two-electron transfer originating from a lone pair of the nucleophile to an antibonding orbital of the metal-oxo species, which is in accordance with the DFT-MD results.

Conclusion and Outlook
This short review on transition metal-based water oxidation focuses on dynamic approaches for a thorough investigation of Comput. 2020, 16, 2436-2449. Copyright 2020 American Chemical Society. [42] water oxidation catalysis by means of DFT-MD and state of the art enhanced sampling methods. In more detail, several studies have highlighted the importance of including environmental effects such as solute-solute interactions and ambient temperature effects when attempting to model properties of proposed intermediates and possible MEPs of different catalysts. To this extent, we have concisely summarized the fundamentals of two selected enhanced sampling methods that can be used to study such systems. We then outlined how different approaches were used to investigate a family of Ru-based WOCs, showcasing potential applications that range from predicting important properties such as pK A values to identifying complex reaction networks and their associated MEPs. The dynamical treatment of the explicit solvation shell thereby enabled us to understand the subtle interplay between the ligand framework and the solvent, which potentially influences the potency of the WOCs at hand. This review also emphasizes the thought process behind our studies and highlights some of their key aspects. In the last section we discussed the O-O bond formation in terms of electron-transfer reactions, where we addressed how the potential failure of DFT can be alleviated.
Over the last two decades we have seen a dramatic increase in computational resources due to the steady development of the required hardware infrastructure. These resources made it possible to study the DFT-MD of reasonably large systems (500-1000 atoms) on a regular basis, circumventing several shortcomings of standardly used static DFT calculations. However, DFT, the workhorse of computational chemists nowadays, possesses limitations when it comes to describe complex multiconfigurational systems. While wave function-based methods can circumvent these limitations, we still face a daily challenge to find the optimal compromise between system size and accuracy. To this end, it is up to the computational chemist to choose an appropriate method and system to tackle the scientific problem at hand.  [29]