Reducing Systematic Uncertainty in Computed Redox Potentials for Aqueous Transition-Metal-Substituted Polyoxotungstates

Polyoxometalates have attracted significant interest owing to their structural diversity, redox stability, and functionality at the nanoscale. In this work, density functional theory calculations have been employed to systematically study the accuracy of various exchange–correlation functionals in reproducing experimental redox potentials, U0Red in [PW11M(H2O)O39]q− M = Mn(III/II), Fe(III/II), Co(III/II), and Ru(III/II). U0Red calculations for [PW11M(H2O)O39]q− were calculated using a conductor-like screening model to neutralize the charge in the cluster. We explicitly located K+ counterions which induced positive shifting of potentials by > 500 mV. This approximation improved the reproduction of redox potentials for Kx[XW11M(H2O)O39]q−x M = Mn(III/II)/Co(III/II). However, uncertainties in U0Red for Kx[PW11M(H2O)O39]q−x M = Fe(III/II)/Ru(III/II) were observed because of the over-stabilization of the ion-pairs. Hybrid functionals exceeding 25% Hartree–Fock exchange are not recommended because of large uncertainties in ΔU0Red attributed to exaggerated proximity of the ion-pairs. Our results emphasize that understanding the nature of the electrode and electrolyte environment is essential to obtain a reasonable agreement between theoretical and experimental results.


■ INTRODUCTION
Polyoxometalates (POMs) are a large group of discrete, polynuclear metal-oxo clusters comprising early transitionmetal (addenda) and oxide atoms. 1−5 Addenda atoms are fully oxidized to d 0 electron configurations capable of forming various topologies employing {MO x } as the principal building block. 1−5 Poly-oxo clusters are formed from the acidification of aqueous molybdate or tungstate oxoanions. 6,7 Partial hydrolysis of polyoxo clusters, achieved through the controlled addition of base produces lacunary clusters. These complexes formally lose one or several M = O vertices and possess reactive cavities with high charge density around the defective region due to the negatively charged oxygen ligands. 8 The defective region can react with transition-metal cations forming a new class of compounds, an example shown are mono-transition-metal-substituted polyoxotungstates, [XW 11 M(L)O 39 ] q− (X = e.g., P(V), Si(IV), L = H 2 O, DMSO, etc.)�see Table 1. 8 Conventionally, hexacoordinate transition-metals are introduced so the sixth coordination site are occupied by solvent ligands from the local environ-ment�see Figure 1. These structures have attracted significant interest as single atom catalysts 9 in fields including water oxidation, 10,11 carbon dioxide reduction, 12,13 and nitrogen activation (Scheme 1). 14 Early computational work on POMs have primarily focused on charged molecules using implicit solvation models. 38 In this regard, Poblet and co-workers reported the relative stability of rotational isomers of Keggin heteropolyanions, α/β-[XM 12 O 40 ] q− . 39 These calculations were performed using the local-density approximation functional coupled with Vosko− Wilk−Nusair parametrization, in the gas phase. 39 The enhanced stability of the α-[XM 12 O 40 ] q− isomer was attributed to the higher energy of the lowest occupied molecular orbital. 39 However, it was shown that four-times-reduced clusters [PW 12 O 40 ] 7− and [SiMo 12 O 40 ] 8− demonstrated enhanced stability (ca. 0.4 eV) toward the β-isomer. 39 Later, Zhang and co-workers emphasized the importance of implicit solvation for reproducing experimental geometries. 40 Several generalized gradient approximation (GGA) functionals were tested, in which Perdew−Burke−Ernzerhof (PBE) and Becke 1988 exchange and Perdew 86 (BP86) methods provided the closest description to the experimental geometries. 40 However, GGA functionals often poorly describe electron delocalization as previously reported by Poblet and co-workers. 41 Their work showed GGA functionals incorrectly localized an additional electron at the belt region in mono-substituted Wells−Dawson [P 2 W 17 MO 62 ] q− (M = V, Mo) anions. 41 On the other hand, employment of hybrid methods, such as B3LYP or M05 correctly localized the electron. 41 However, only B3LYP [20% Hartree−Fock (HF) exchange] gave the correct ordering and relative reduction energies with respect to experimental measurements. 41 Quantum chemical calculations reproducing redox potentials in POMs have been reported without accounting for electrolyte environment and counterions, using implicit solvation models.
Aparicio and co-workers computed the tungsten redox waves for mono-substituted Keggin, [XMW 11 O 40 ] q− (M = W, Mo, V, Nb, and Ti) using a conductor−like screening mode. 42 This simplification produced large uncertainties in absolute reduction potential. However, shifts in potential with respect to the Keggin [XM 12 O 40 ] q− anion were reproduced with corresponding experimental data. 42 Recently, Rosch and co-workers reported redox potentials for the tri-Mn-substituted Keggin in which counterions were explicitly located onto the surface of POM to neutralize the system. 43 The authors employed several exchange−correlation (x−c) functionals with the increase in contributions of HF exchange (0% PBE, 10% TPSSh, 20% B3LYP, 25% PBE0) which produced redox potentials that increase (in that order) with the exact exchange. 43 The closest agreement with experimental literature was shown for the GGA−PBE functional attributed to fortuitous error cancellations. By contrast, hybrid exchange−functionals overestimated experimental potentials by 0.6−1.0 V. 43 This work was extended to include explicit treatment of water molecules derived from molecular dynamics simulations. 44 The inclusion of explicit solvation significantly improved computed redox potentials; however, its contributions were out-weighted by the choice of exchange−correlation functional. 44 In this work, optimal reproduction of experimental literature were obtained using the B3LYP method. 44 Recently, Falbo and Penfold reported the influence of the self-interaction error could be reduced by neutralizing the charged species by incorporating counterions into the system. 45 The authors reported excellent agreement, within 0.1 V of the experimental literature for the first redox wave in Na 3 [SiW 12 O 40 ]. 45 In this work, we have performed systematic density functional theory (DFT) calculations to study the redox properties of [PW 11 M(H 2 O)O 39 ] q− M = Mn(III/II), Fe(III/II), Co(III/II), and Ru(III/II) and their potassium salts. For achieving adequate molecular geometries, we present a structural benchmark for the cobalt(II)-substituted Keggin against the crystallographic structure. 46 Thereafter, we explore the challenges in computing redox potentials and provide an insight into the geometric and electronic factors controlling it.

■ COMPUTATIONAL DETAILS
All computational results were obtained using the ARCHIE− WeSt high-performance computer based at the University of Strathclyde. DFT calculations were performed using the Amsterdam Modelling Suite (AMS 2020.1) package. 47 In this work, several classes of exchange−correlation (x−c) functionals were employed, which include (i) generalized gradient approximation (GGA); (ii) hybrid; and (iii) range-separated hybrid functionals. GGA functionals considered were as follows: (i) PBE; 48 (ii) Perdew−Wang (PW91); 49 and (iii) Becke 1988 exchange and Perdew 86 (BP86). 50,51 The hybrid x−c functionals considered were as follows: (i) Becke, 3-parameter, Lee−Yang−Parr (B3LYP*, 52 B3LYP 53 ); (ii) PBE0; 54 and (iii) Becke's half-and-half (BH&H). 55 Hybrid functionals were selected on their contributions of HF exchange (15% B3LYP*, 20% B3LYP, 25% PBE0, and 50% BH&H). The ωB97X method was selected as the range-separated hybrid functional. 56 We employed Slater basis sets comprising the following: (i) triple-ζ polarization (TZP); (ii) triple-ζ plus polarization (TZ2P); and (iii) quadruple-ζ plus polarization (QZ4P). 57,58 Relativistic corrections were included by means of the zeroth order regular approximation formalism. 59 The effects of aqueous solvent were approximated by using the conductorlike screening model, as implemented by AMS. 60 For open shell molecules, unrestricted Kohn−Sham (UKS) theory was implemented, while restricted Kohn−Sham (RKS) theory was employed for closed shell systems. All harmonic vibrational frequencies were calculated using PBE coupled with the TZP basis set. The calculation of Gibbs free energies for hybridoptimized systems were corrected by using the zero-point energies and entropic components obtained from GGA− vibrational frequencies�see eq 1.0.  To evaluate the discrepancy of the calculated versus crystallographic geometries, we employed mean absolute error (MAE), mean signed error (MSE), and standard deviation (STD) calculated using where d calc and d exp are the calculated and experimental bond distances, respectively.   46 However, our models have inverted this observation, in which the magnitude of O c −Co exceeds O t −Co. Across all methods, the equatorial distances were accurately described with computed errors rarely exceeding 0.02 Å. The equatorial parameters were insensitive across all applied functionals and basis sets, reporting a range of 0.059 and 0.085 Å for O a1 −Co and O b2 −Co, respectively.
The electronic structure of a fully oxidized POM consists of two identifiable bands: (i) an oxo-band comprising occupied oxo-ligands; (ii) and a metallic band encompassing unoccupied addenda orbitals. Figure 3 23,24 and Ru(III/II). 24,26,27 Figure 3 and Table-S1 39 ] q− were calculated without neutralizing their charge, and therefore, large discrepancies associated with the self-interaction error are expected for these calculations. Previous work showed explicitly locating counterions on the surface of the POM, rendering the system charge neutral can reduce this systematic error. 44

Computation of Redox Potentials-Charge Neutral
Model. On the surface of the molecule, the counterion can interact with 4-fold (pocket A−F) or 3-fold pockets�see Figure  6. The 4-fold pockets feature four oxygen atoms capable of interacting with counterions, while 3-fold pockets are comprised three adjacent oxygen atoms which assume a triangular shape. In this work, we concentrated on counterion interactions with the 4-fold pockets, which bind significantly stronger than the 3-fold pockets due to improved Coulomb interactions.

Inorganic Chemistry pubs.acs.org/IC Article
To differentiate between cation arrangements, isomers were identified by their unoccupied pockets. For example, isomer-D corresponds to systems with an absent cation−oxygen interaction at pocket D�see Figure S3. The relative stability of all isomers was evaluated by the difference in electronic energy, E Rel .�see isomers were almost negligible at 1.11 kcal mol −1 . The discrepancy between pocket A and B was attributed to reduced electrostatic attractive forces because of shielding by the aqua protons, correlating with previous work. 43  Hence, it is reasonable to assume all arrangements will co-exist in solution. However, for our redox calculations, we focused on the two lowest energy isomers: A and A,D.
Redox potentials, U 0 Red versus SHE, for Mn(III/II), Fe(III/ II), Co(III/II), and Ru(III/II) couples present in K x [XW 11    Error computed using BH&H ranged from 0.43 to 1.50 V. Hybrid functionals exceeding 25% HF exchange did not provide any improvement for reducing ΔU 0 Error . As is evident, increasing HF exchange positively shifted U 0 Red which was attributed to the over-stabilization of the ionpairs. Hybrid functionals exceeding 25% HF are not recommended because of large CPU times coupled with significant overestimations to ΔU 0 Red . U 0 Red became positively shifted with the increase in contributions of HF exchange. To rationalize this observation, we plotted counterion−bridging oxygen (between pocket D) as a function of HF exchange, see Figure 10. Increasing contributions of HF exchange produced shorter O b −K distances leading to the over-stabilization of the close contact ion-pairs. The significance of these O b −K distances is emphasized by the computed range (PBE to BH&H) in O b −K distances in K 5 39 ] which reflected in the range of 1.29 V for U 0 Red . Obtaining accurate redox potentials using implicit solvation models remains a challenge. Our results emphasize that understanding the nature of the electrode and electrolyte environment may be essential to obtaining reasonable agreement between theoretical and experimental results. Further improvements to this work by accurately modeling counterion− bridging oxygen distance will enable more accurate computation of U 0 Red . The present results emphasis the current approach requires significant improvement to achieve more reliable modeling with respect to experimental work.