A Force Field for a Manganese-Vanadium Water Oxidation Catalyst: Redox Potentials in Solution as Showcase

: We present a molecular mechanics force ﬁeld in AMBER format for the mixed-valence manganese vanadium oxide cluster [Mn 4 V 4 O 17 (OAc) 3 ] 3 − —a synthetic analogue of the oxygen-evolving complex that catalyzes the water oxidation reaction in photosystem II—with parameter sets for two different oxidation states. Most force ﬁeld parameters involving metal atoms have been newly parametrized and the harmonic terms reﬁned using hybrid quantum mechanics/molecular mechanics reference simulations, although some parameters were adapted from pre-existing force ﬁelds of vanadate cages and manganese oxo dimers. The characteristic Jahn–Teller distortions of d 4 Mn III ions in octahedral environments are recovered by the force ﬁeld. As an application, the developed parameters have been used to calculate the redox potential of the [Mn III Mn IV3 ] (cid:10) [Mn IV4 ]+e − half-reaction in acetonitrile by means of Marcus theory.


Introduction
There is a large consensus among scientists that global warming threatens life on Earth as we know it [1]. To lower CO 2 emissions that contribute to global warming, fossil fuels need to be replaced by sustainable energy sources such as sunlight [2]. Artificial photosynthesis-photochemically driven water splitting [3]-is emerging as a promising and clean technology to convert sunlight into a storable energy form. Devices for this task consist of several components, where the water oxidation catalyst is among the most critical ones [4][5][6][7][8][9].
In natural photosynthesis, the oxidation of water oxygen is enabled by a metalloenzyme [8], called the oxygen-evolving complex (OEC). The active center of the OEC contains an oxo-manganese, cubane-like structure with a Ca 2+ ion occupying one of the vertices. Due to the extremely high efficiency displayed by this catalyst (with a turnover number of about 600,000 and a turn-over frequency of 30-90 s −1 according to the literature [10,11]), there have been major efforts in attempting to understand the mechanism of the overall water oxidation catalytic cycle, including the nature of the intermediates that participate in that mechanism [12][13][14][15][16].
The outstanding efficiency of the natural OEC makes it an appealing model for the development of analogous artificial water oxidation catalysts (WOCs). Based on this idea, several synthetic WOCs based on metal oxide clusters-polyoxometalates (POMs)-have shown promising catalytic activities towards the water oxidation reaction [5][6][7][8]10]. In these POM-based WOCs, typical transition metal complexes are either rare transition metals like Ru or Ir [5], or earth-abundant metals like W, Mo, Ni, Co, Fe, or Mn [7].

Theory
This section documents the main working equations used in the parametrization of the FF for the MnV WOC and for the computation of the reduction potential.

Force Field Parameters
The general FF form implemented in the AMBER [24] package consists of a sum of harmonic potentials for bond stretching and angle bending, plus a truncated Fourier series for describing dihedral torsions, and the nonbonded Lennard-Jones and pairwise electrostatic potential functions, as described by the following equation [24]: This expression is explained in more detail in Section S1 in the Supporting Information (SI).
In the present work, the initial guesses for the harmonic bond and angle parameters (r i,eq , K r,i , α i,eq , K α,i ) were obtained by using data from quantum mechanics/molecular mechanics (QM/MM) reference trajectories (see Section 3.2), and fitting each of the bond and angle distributions obtained to a Boltzmann distribution of the form: where N is a normalization constant, k B is the Boltzmann constant, and T is the temperature.
The energy E i can be represented with the harmonic oscillator expression Here, X i,eq indicates either the bond or the angle equilibrium value (r i,eq , α i,eq ), depending on the considered parameter. Within this framework, each bond (angle) distribution is, in practice, fitted to a Gaussian curve, so that the equilibrium bond (angle) would correspond to the average of the distribution and the force constant K i would be associated with its width (i.e., the standard deviation σ) as We note that this approach neglects a factor arising from the Jacobian of the Cartesianto-internal coordinate transformation, see, for example, Ref. [25]. However, as we use Equations (2)-(4) only to obtain the initial parameter set (before refinement), omission of this factor does not affect the accuracy of our final parameters.

Redox Potentials
In the present work, we focus on the [Mn III Mn IV 3 ]/[Mn IV 4 ] redox couple, whose oxidation process is given by which can be schematically represented as in which R and O represent the reduced ([Mn III Mn IV 3 ]) and the oxidized ([Mn IV 4 ]) states, respectively. In general, the redox potential of the reaction shown in Equation (6) is determined from the Gibbs free energy of reaction ∆G R→O as follows: where ∆G ref is the Gibbs free energy difference for the redox reaction in the reference electrode, n is the number of exchanged electrons, and F is the Faraday constant. The approach adopted in the present work for calculating ∆G R→O is based on the methodology first proposed by Warshel [26,27] and applied to one-electron half-reactions of complex systems in the last few years (e.g., see Refs. [28][29][30]). The starting point is the following equation: in which ∆A R→O is the Helmholtz free energy of oxidation, p is the pressure, and V is the volume of the system. In the case in which the oxidation reaction is followed by small structural changes, the volume difference can be neglected, and therefore the Gibbs free energy ∆G R→O (obtained at constant pressure) can be assumed to be equal to the Helmholtz free energy ∆A R→O (obtained at constant volume). Although ∆A R→O , in principle, depends on the canonical partition functions for the oxidized and the reduced states, it can be shown that ∆A R→O can be calculated from the distribution of vertical energy gaps between the two states [27]: Here, {R N } is a geometry sampled from a given trajectory. Thus, depending on whether {R N } corresponds to a geometry of the reduced or of the oxidized species, ∆E R→O ({R N }) indicates the ionization energy (IE) or the electron affinity (EA) correspondingly. With this expression for ∆E R→O ({R N }), ∆A R→O is given by [31]: where R indicates an average over all geometries sampled from the trajectory of the reduced state, and O is the analogue for the oxidized state. This expression can be further simplified by applying Marcus theory [32][33][34], in which it is assumed that the solvent responds linearly to a change in the solute. As a consequence, the distributions of the vertical energy gaps ∆E R→O ({R N }) for both the reduced and the oxidized species-that is, the distributions of the IEs and the EAs, respectively-will be Gaussian functions, and their standard deviations will be the same. Hence, the Helmholtz free energy difference expression (Equation (10)) can be simplified as follows: Equation (11) is the working equation for the determination of the Gibbs free energy of oxidation in the present work. ∆E R→O R and ∆E R→O O are the average IE and the average EA, respectively. Within the linear response regime, the reorganization energy λ from Marcus theory can also be computed from the average IE and EA [28][29][30]: Additional estimates of the reorganization energy for R and O can be obtained directly from the standard deviations of the IE and EA distributions: where σ R or σ O are the standard deviations of the distributions of vertical energy gaps for the R trajectory (IE distribution) or for the O trajectory (EA distribution). If the three estimates for the reorganization energy (λ, λ R , λ O ) do not agree, this indicates a deviation from the Marcus regime. Using the "Q-model" developed by Matyushov and Voth [35,36], one can compute an error estimate for the free energy [30]:

Computational Details
In this section, we describe the details of the different steps performed. The first step consists of ab initio calculations to optimize the equilibrium geometries and extract partial charges. The second step is the construction of the system from solute, solvent, and counter ions, and to compute two hybrid QM/MM molecular dynamics (MD) trajectories as references to fit the FF parameters. The parameters were manually improved and iteratively based on the comparison between the QM/MM MD trajectory and MM MD trajectories generated with the current FF parameters. With the final set of parameters, the reduction potential was computed.

Reference Ab Initio Calculations
The geometries of [Mn III Mn IV 3 ] and [Mn IV 4 ] were optimized at the U-BP86 [37,38] level of theory (a high-spin configuration was adopted throughout), using the ZORA-def2-TZVP basis set for the description of the Mn, V, O atoms, and the ZORA-def2-SVP basis set for the C and H atoms [39,40]. The ZORA relativistic Hamiltonian [40][41][42] was used to account for scalar relativistic effects. Empirical dispersion correction effects were introduced through Grimme's D3 correction [43]. Tight convergence criteria were used throughout the geometry optimizations. The Mulliken charges of the final converged calculation were used as the partial charges within the FF. These calculations were performed with the ORCA 4.0.1 [44,45] software.
The combination of U-BP86, ZORA, and D3 dispersion correction was previously shown to describe geometries of first-row transition metal complexes very well [14,46,47]. Furthermore, as a GGA functional, U-BP86 allows for very efficient calculations that are beneficial for the expensive QM/MM reference trajectories. We note that this level of theory is, in principle, not consistent with the electronic structure used for parametrizing the generalized AMBER force field (GAFF) [48] that we use for describing the solvent. However, the GAFF methodology-which uses Hartree-Fock and Møller-Plesset second-order perturbation theory-was never designed for complicated inorganic compounds like POMs, and in fact, preliminary Hartree-Fock calculations did not even converge for the MnV WOC. Hence, we believe that a description of this molecule through DFT will lead to more accurate FF parameters and also to more realistic solute-solvent interactions than an approach using Hartree-Fock or Møller-Plesset perturbation theory.

QM/MM MD Reference Simulations
In order to refine the FF parameters for [Mn III Mn IV 3 ] and [Mn IV 4 ], we first produced two QM/MM MD trajectories to serve as a reference. These simulations were performed by using the AMBER 17 [49] package for describing the nuclear motion, along with its interface with the Terachem [50,51] package for the electronic structure calculations. The total system consisted of the MnV WOC, one (two) Na + counterion(s) for the oxidized (reduced) species, and 2300 acetonitrile molecules in a periodic cubic box with a side length of about 58 Å. This solvent was chosen in order to mimic the experimental conditions [17] that employ dry acetonitrile as a solvent for the electrochemical measurements. Before the QM/MM simulations, the system was first minimized, heated to 300 K (NVT ensemble, Langevin thermostat), and equilibrated to a pressure of 1 bar (NPT ensemble, Berendsen barostat) over 120 ps using MM MD with the initial set of FF parameters (see below) and by imposing a force constant of 200 kcal mol −1 Å −2 on the MnV WOC. This equilibration was followed by a 5 ps QM/MM unconstrained equilibration and a production phase of another 5 ps, using the same conditions. From the last 5 ps, we extracted 100 equidistant snapshots to generate the distributions of bond lengths and angles for use with Equation (4) and as a reference for the refinement of the harmonic terms. A time step of 0.2 fs was used for classical MD runs, whereas a 0.5 fs time step was used for the QM/MM calculations.
In the QM/MM MD calculations, the QM region contained the [Mn III Mn IV 3 ] or [Mn IV 4 ] molecule (46 atoms) and was treated with the U-BP86 [37,38] level of theory. The Mn and V atoms were described by the LANL2DZ effective core potential basis set [52][53][54], whereas the C, H, and O atoms were described by the def2-SVP [39] basis set. The acetonitrile molecules and the counterions were described with GAFF [48] in the same way as in the classical MM MD simulations, see below. It was previously shown that GAFF describes this particular solvent very well [55,56]. In order to account for the electrostatic interactions between the QM nuclei and the MM point charges that enter the Hamiltonian, a cutoff of 10 Å was applied.

MM MD Simulations
The same setup of the system as the one described above for the QM/MM MD simulations has been used for the classical MM MD simulations. These simulations have been performed using the GPU accelerated version of the AMBER 17 [24] package, again using the GAFF force field [48] for the acetonitrile solvent molecules and the Na + counterions, and the newly parametrized FF for [Mn III Mn IV 3 ] or [Mn IV 4 ]. Minimization, heating, and equilibration were performed as described above. A 500 ps unconstrained equilibration was performed prior to the 10 ns production phase, from which 100 snapshots were sampled (every 100 ps) to compute bond lengths and angles, which were then compared to the QM/MM distributions for refining the FF parameters. A time step of 0.2 fs was used throughout. The 100 snapshots of the final MM MD trajectories were also used to compute Mulliken spin populations (to analyze the stability of the electronic wave function), using the same level of theory as used in the QM/MM MD trajectory.

Parameter Setup
The FF presented in this paper reuses a few atom types and parameters from GAFF [48] (although we do not follow the GAFF parametrization methodology because it is not designed for describing POMs), whereas the parameters involving Mn, V, and O atoms were specifically adapted to the MnV WOC. Here, we focus our attention on finding reasonable bond and angle parameters for the [Mn III Mn IV 3 ] and [Mn IV 4 ] molecules. As in these complexes, the bonds and angles involving the Mn atoms are strongly dependent on the electronic structure of the molecule (through the Jahn-Teller effect), we estimate that many of the provided parameters require careful reconsideration when adapting our FF parameters for similar POM molecules.
The following choices for the FF parameters were made in this work. (i) The partial charges q i for the MnV WOC were obtained as Mulliken charges from the ab initio calculations (Section 3.1). (ii) The Lennard-Jones parameters were taken from AMBER for Mn, O, C, H, and Na + [57,58], whereas the parameters for V were adapted from Ref. [59]. (iii) Dihedral terms do not play a large role due to the rigidness of the polyoxometalate structure and were adapted from the literature [59,60] or taken from GAFF for the acetate ligands.
Initial guesses for the bond and angle parameters (iv) of the molecule were either taken from the literature [59,60] or from the bond length/angle averages of the QM/MM MD trajectory. The force constants for bonds were computed from the QM/MM bond distributions and Equation (4). The angle force constants were kept at the initial values. With the full set of parameters (i-iv), an MM MD run was performed as described in Section 3.3, and the bond and angle distributions computed. From a comparison of the QM/MM MD and MM MD distributions, the bond and angle parameters were refined, and a new MM MD simulation was carried out. This procedure was repeated iteratively until the root-mean-square deviations (RMSDs) of the average bond lengths and angles (between MM MD and QM/MM MD) were below 0.02 Å or 2 degrees, respectively.
The final set of parameters is presented in Section S1 in the SI, and is also given in the Supplementary Materials in AMBER format. The absolute deviations of average bond lengths and angles between QM/MM MD trajectories and MM MD trajectories are presented in Section S2 in the SI.
We note that besides manual refinement of FF parameters, as carried out here, there are also a large number of automatic procedures available in the literature, for example, the Seminario method [61]. Although the prospect of automatic parametrization is promising, preliminary tests showed that the Seminario method cannot be easily applied to the MnV WOC with its complicated and uncommon connectivity. The FF parameters that we obtained in this way performed significantly worse than the manually refined parameters, which is why we focus only on the latter in this work. Other possible alternatives to obtaining FFs for challenging systems include QM-derived FFs [62,63] and machine-learning-based potentials [64,65].

Redox Potential Calculations
To compute the redox potential of the [Mn III Mn IV 3 ]/[Mn IV 4 ] pair, we performed two more MM MD simulations. We used the same setup and preparation steps as described above, but the production run was extended to 100 ns (0.2 fs time step). The MnV complex was described by the final versions of the newly developed FF. From each of the two trajectories, we extracted 2000 snapshots (at 50 ps intervals).
At each of the 4000 snapshots, we carried out single point calculations for the [Mn IV 4 ] wave function (charge −1 and multiplicity 13) and the [Mn III Mn IV 3 ] wave function (charge −2 and multiplicity 14). These calculations were both carried out at the U-BP86 [37,38] and U-B3LYP [66] levels of theory, employing the ZORA-SVP basis set to reduce the high computational cost. As before, we applied the ZORA relativistic Hamiltonian and the D3 correction. The surrounding point charges were included in the calculation through the AMBER-ORCA QM/MM interface [67]. Here, no periodic boundary conditions were applied, and all MM point charges in the primary simulation cell (13801 or 13802 point charges per calculation) were included in all calculations. Prior to the single point calculations, the solute was imaged to the center of the simulation cell. For the interested reader, the strong influence of the point charge placement on the obtained energy gaps is documented in Section S3 in the SI.
Finally  (12) and (13). The final value of the reduction potential was computed using Equation (7), where ∆G ref was 4.998 eV, which is the absolute electrode potential of the ferrocenium/ferrocene (Fc + /Fc) reference electrode [68]. Error estimates were obtained from Equation (14). The two QM/MM MD trajectories that were simulated were not employed for the reduction potential calculations. As documented in Section S4 in the SI, the short simulation time of the QM/MM MD trajectories precludes proper sampling of the solvent around the MnV WOC. Whereas the 100 ns MM MD trajectories can uniformly sample the first solvation shell, in the 5 ps QM/MM MD trajectories, the solvent is restricted to configurations that are too close to the final geometry from the equilibration run. Preliminary calculations showed that this essentially "frozen" solvent leads to a strong bias in the computed energy gaps, and therefore to a very inaccurate redox potential.

Results
In the following, we discuss the implications of the catalyst structure on the atom type labeling, the quality of the force field with respect to accurate bond lengths and angles, and the computed redox potentials.

Structure of the Catalyst and Atom Type Labeling
The full three-dimensional structure of the MnV WOC is shown in Figure 1a  The entire molecule consists of 46 atoms, but like in many other FFs, not every atom necessarily requires its own independent set of FF parameters. Hence, here we use the concept of "atom types" like in GAFF [48] to define our FF parameters, where atoms with the same atom type label (e.g., "MA" or "c3") use the same parameters. The assignment of these atom type labels is shown in Figure 1c-e. Here, to avoid occlusion, we only show part of the connectivity network in each of the three panels. The panels can be obtained by a rotation of the entire molecule by 120 degrees (or 240 degrees) before the occluded atoms are deleted. Note that some of these atom types were directly reused from GAFF (atom types c, c3, hc, and o in the acetate ligands) and those atom type labels are written in lower case, as is customary in GAFF. Conversely, all newly defined atom type labels are upper case to not clash with the ones already present in GAFF.
The [Mn IV 4 ] structure presents an idealized C 3v symmetry, where the rotation axis passes through the apical Mn, the opposite oxygen atom, and the central V=O group. This ideal symmetry of the reactant is, however, dependent on the local oxidation states of the manganese atoms. In particular, Mn III centers have a d 4 configuration that leads to a characteristic distorted octahedral coordination induced by the Jahn-Teller effect [14]. This will lead to strongly differing bond lengths for Mn III and Mn IV , thus lowering the symmetry of the system. To support lower symmetry than C 3v in our topology, we have assigned different atom types to the Mn and O atoms in the cubane, and to most atoms in the vanadate. As the acetate geometry is virtually unaffected by the oxidation state of the cubane, different atom types for the acetates were not necessary and we directly reused the lower case GAFF atom types c, c3, hc, and o.
The Mn atoms are labeled as MA, MB, MC, and MD, where MB is the apical Mn atom bound to all three acetate ligands, and MA, MC, and MD are ideally symmetry equivalent. The cubane oxygen atoms are labelled as OJ, OK, OL, and OM, where OJ is opposite to the apical MB and OK, OL, and OM are symmetry equivalent. The apical V atom is VA, the other ones are VB and VC. The vanadate oxygens were labeled OA, OC, and OQ for the V=O groups; OB and OR for the V-O-V bridges; and OD, OE, OF, OG, OH, OI for the V-O-Mn bridges. For the acetates, we have kept the GAFF nomenclature, except that we assigned new labels (ON, OO, OP) to the oxygen atoms bound to MB, in order to enable unequal bond parameters with that atom (i.e., in order to avoid three MB-o bonds that could not be parametrized independently).

Force Field Parameters and Geometries
In order to check the accuracy of the final FF parameter set, we present the averages and standard deviations of all bond lengths and angles in (largest deviation 7.7 degrees). This is sufficiently good for a qualitative (and possibly quantitative) description of the dynamics of the MnV WOC, as we will scrutinize below.
Notably, most of the average bond lengths calculated from both the MM MD and the QM/MM MD trajectories are in good agreement with the experimental values of the [Mn III 2 Mn IV 2 ] pristine catalyst obtained from X-ray crystal diffraction [17], even though they correspond to different oxidation states. This is shown in Figure 2b Panels (a,b) show the averages and standard deviations of the bond lengths, whereas panels (c,d) show the corresponding data for the bond angles. In (b,d), rectangles show the bond lengths from X-ray crystallography for comparison [17]. Bond lengths discussed in the text are marked with arrows in (b).
In Figure 3 we present the Mulliken spin populations of the [Mn III Mn IV 3 ] complex at the 100 snapshots along the reference QM/MM MD trajectory and the MM MD trajectory with the final parameter set. As can be seen, at our level of theory, the MnV WOC shows a valence localization of the oxidation states, where the atoms MA, MB, and MD are in oxidation state +IV (indicated by about 2.9 unpaired electrons) and atom MC is in oxidation state +III (about 3.7 unpaired electrons). The different oxidation states are the reason for the different bond lengths shown in Figure 2b; the MC-o and MC-OJ bonds (indicated by the left-most arrows) are strongly elongated, whereas the equivalent bonds of MA and MD are not. Interestingly, the oxidation states remain localized for the entire length of the QM/MM MD trajectory (including the QM/MM part of the equilibration) over 10 ps. For the much longer MM MD trajectories, the oxidation states also remain localized, although it can be expected that this is due to the fact that the FF parameters enforce the Jahn-Teller axis to remain at OJ-MC-o.
The finding of localized electronic states is in agreement with results on comparable manganese oxygen clusters like the natural OEC [14,69]. The reason why the crystal structure indicates delocalization of the oxidation state (evidenced by the equal OJ-Mn-o bond lengths) could be that there is a small activation barrier that needs to be overcome to move an electron from Mn III to Mn IV . Thus, on microscopic time scales, there is valence localization, but over macroscopic time scales, the three equivalent Mn atoms are, on average, in oxidation state (3 + 4 + 4)/3 = 3.67 in [Mn III Mn IV 3 ], with the apical Mn atom being formally in oxidation state 4. This agrees reasonably well with the results of Schwarz et al. [17] using the bond valence sum method [70,71], who assigned (for [Mn III 2 Mn IV 2 ]) a value of 3.4 to the MA, MC, and MD atoms and 4.3 to the MB (apical) atom.  The angle parameters were tuned in a similar manner as the bond parameters, but only the equilibrium angles were refined with the average values of the QM/MM MD reference trajectories. As presented in Figure 2c,d, there is very good agreement between the QM/MM MD and the MM MD trajectory, showing differences of at most 7.8 degrees, a value which is within the order of magnitude of the standard deviations observed. The largest deviations from the X-ray diffraction values (rectangles in Figure 2d) can again be explained by the valence localization of the Mn III oxidation state and the ensuing Jahn-Teller distortions. The good agreement of the angles after adaptation of the equilibrium angle parameters made it unnecessary to further adjust the angle force constants, which were therefore kept as used the initial set of parameters (Table S2 in the SI).
As shown in Figure 2, in most cases the standard deviation (the width of the distribution of bond lengths or angles) is smaller than the QM/MM MD counterpart, despite the refinement procedure described above. Although the force constants could in principle be manually modified in an attempt to match the broadness of the bond length distribution obtained, this is not trivial, as smaller force constants will make the bonds weaker and thus will likely lead to larger deviations of the averages. To match both averages and standard deviations would likely require some numerical or analytic procedure that also compensates for the intramolecular Coulomb interactions, like the force matching procedure [72,73]. Additionally, the narrower distributions in the MM MD trajectory might be partially due to the use of harmonic spring forces in the FF, which is suitable for small displacements but fails to describe the anharmonic potential energy surface for larger displacements [74].

Energy Distributions and Redox Potentials
The calculation of the redox potential for the redox couple [Mn III Mn IV 3 ]/[Mn IV 4 ] is based on two 100 ns MM MD trajectories (one for the reduced and one for the oxidized species), from which 2000 snapshots each were sampled and single point SCF calculations were performed. Prior to the single point calculation, the solute was positioned at the center of the simulation cell and then all point charges from the primary cell were considered in the computation. The vertical ionization energies and electron affinities were computed as energy differences between the SCF energies of the two single point calculations at the same geometry.
The distributions of vertical energy gaps obtained from the single point calculations are shown in Figure 4, whereas the derived results are collected in Table 1. The absolute redox potential ∆A R→O computed from Equation (11) was +5.59 eV for BP86 and +6.31 eV for B3LYP. The reorganization energy can be computed in three different ways, either from the half difference of the average energy gaps (λ), or from the standard deviations of the energy gaps (λ O and λ R ). The degree of agreement between these three values offers a way to judge whether the reaction follows the linear response Marcus regime and, hence, whether the equations given in Section 2 are applicable. As the different reorganization energy estimates differ among themselves, we used the so-called Q-model approach [35,36] to verify our results. The error estimates from this model, Equation (14), are satisfyingly small, being only 0.015 eV for BP86 and 0.025 eV for B3LYP. The final reduction potential of the [Mn III Mn IV 3 ] [Mn IV 4 ]+e − half reaction is then computed from Equation (7), using the value of +4.998 eV as the absolute redox potential of the Fc + /Fc reference electrode [68] that was employed by Schwarz et al. [17]. The corresponding uncertainty was estimated as Thus, the reduction potential of the [Mn III Mn IV 3 ] [Mn IV 4 ]+e − half-reaction is found to be +0.59 ± 0.28 V (BP86) or +1.31 ± 0.30 V (B3LYP). These values can be compared to the experimental value of +1.1 eV determined by cyclic voltammetry with the same reference electrode [17]. Here, multiple error sources might explain the differences between experimental and computed values. First, the neglected p∆V R→O term in Equation (8), is probably insignificant, as the box volumes of the two MM trajectories are identical to within 0.03 standard deviations. However, there might be small contributions to the uncertainty from entropic contributions due to the triply degenerate Jahn-Teller configurations that are not taken into account in our FF and from spin coupling effects.  (11) 5.59 6.31 λ from Equation (12) 1.78 2.02 λ O from Equation (13) 1.26 1.33 λ R from Equation (13) 1.37 1.52 δ∆A corr from Equation (14) −0.015 −0.025 E R→O from Equations (7) and (8) +0.59 ± 0.28 +1.31 ± 0.30 Second, formally the probability distributions in Figure 4 are not entirely correct, as the MM Hamiltonian used for sampling and the QM/MM Hamiltonian used to calculate the energy gaps are not the same, and their probability distributions will differ to some extent [75][76][77]. To check to which extent this can introduce a bias in our results, we investigated how much the MM and QM/MM potential energies deviate from each other and and whether the vertical energy gaps depend on this deviation. As we show in Section S5 in the SI, the MM and QM/MM potential energies show deviations of a few kcal/mol. However, the vertical energy gap distributions do not depend on the potential energy deviations, showing that sampling with the MM Hamiltonians most likely did not introduce a significant bias into the distributions of the vertical energy gaps. Hence, for the sake of statistical convergence of the histograms in Figure 4 we omit unbiasing here, as was also done occasionally in the past when MM-based sampling and QM/MM computation of the energy gaps were combined to compute redox potentials [78]. Further improvement of the sampling could be achieved in future work by relaxing each MM snapshot with short QM/MM trajectories, a procedure that is sometimes called "switching" [77]. This might especially help in alleviating the fact that the harmonic terms in our optimized FF are systematically too stiff.
Third, although the accuracy of the FF might introduce some bias, based on the comparison of our results with two different density functionals-the GGA functional BP86 and the hybrid functional B3LYP-it appears that the electronic structure level of theory has an even stronger influence on the final result. Here, apparently the exact exchange in the B3LYP functional reduces the self-interaction error and thus stabilizes the reduced species relative to the oxidized species, which increases the energy gaps and thus shifts the reduction potential to more positive values. From the comparison with the experimental value, B3LYP provides a more accurate estimate of the reduction potential (BP86 is off by about 0.4 eV, B3LYP by about 0.2 eV). However, it seems that B3LYP overshoots these electronic effects. Hence, some more elaborate density functional, and probably a larger electronic basis set, might be necessary to obtain the right balance between electronic correlation and exchange, although such investigation of electronic structure effects goes beyond the scope of the present work.
Finally, a comment on the importance of proper solvent sampling is on order here. Besides the MM MD trajectories, in preliminary calculations we also attempted to compute the reduction potential from the snapshots of the QM/MM MD trajectories. Using the BP86 functional, the final reduction potential based on a subset of QM/MM MD snapshots was found to be −1.10 ± 0.42 V (against the Fc + /Fc electrode), far away from the experimental value of +1.1 V and the values based on MM MD (+0.59 V or +1.31 V, see Table 1). As shown in Section S4 in the SI, one of the most important differences between the QM/MM MD and MM MD simulations is the accomplished solvent sampling, which is most likely the reason for the large differences in the resulting reduction potentials. This illustrates that proper solvent sampling-requiring long time scales-is mandatory for redox potential computations like the ones presented here. Naturally, such long time scales and uniform solvent sampling can be achieved using MM FFs, showing that an accurate FF is of large importance for future atomistic simulations of the embedding of the MnV WOC in block copolymer matrices.

Conclusions
We have developed a force field in AMBER format for a manganese vanadium oxide synthetic analogue of the oxygen evolving complex, which catalyzes the water splitting reaction in nature. In particular, we reported on two sets of force field parameters for the two oxidation states involved in the [Mn III Mn IV 3 ] [Mn IV 4 ]+e − oxidation reaction that precedes the water oxidation reaction [17] in this catalyst.
Initial force field parameters were obtained by combining pre-existing force fields of polyoxovanadate cages [59], high-valent manganese dimers [60], and ab initio calculations. This force field was refined by tuning the equilibrium values and the force constants of the bond and angle potentials, using the distributions of bonds and angles obtained from reference QM/MM MD trajectories. It was observed that the refined force field parameters properly describe the axial Jahn-Teller distortion characteristic of a d 4 Mn III center in an octahedral ligand field. Furthermore, we found that along the dynamics, according to the Mulliken spin populations and the bond lengths of the bonds involved in that geometrical distortion, no valence change or electron transfer occurred-hinting at valence localization in this complex. In general, good agreement of bond lengths and angles between the MM MD, QM/MM MD, and (with exceptions) the X-ray structure [17] was found.
The redox potential of the [Mn III Mn IV 3 ] [Mn IV 4 ]+e − oxidation reaction was calculated by applying a protocol originally proposed by Warshel and coworkers [26,27], in which the Gibbs free energy of oxidation was approximated as the Helmholtz free energy (Equation (11)), calculated as the average between the average ionization energy and the average electron affinity. Using two different density functionals, we computed the reduction potential to be +0.59 ± 0.28 V (BP86) or +1.31 ± 0.30 V (B3LYP). Although these results are in reasonable agreement with the experimental value of +1.1 eV [17], it appears that the chosen electronic structure level of theory has a strong influence on the results. Importantly, we also found that qualitative agreement with experimental values requires extensive sampling of the solvent shell around the MnV WOC, as it can currently only be achieved with MM force fields but not with QM/MM methods. We thus expect the presented force field parameters to enable reasonable atomistic simulations of challenging polyoxometalate complexes.

Data Availability Statement:
The data that support the findings of this study are available from the corresponding author upon reasonable request.