Calibration of Computational M¨ossbauer Spectroscopy to Unravel Active Sites in FeNC-Catalysts for the Oxygen Reduction Reaction

,


Introduction
2][3][4][5][6][7][8][9] The experimental study of iron in such systems is often complemented with a computational analysis, chiefly via density functional theory (DFT).As a result, the oxidation and spin states of the iron ion, its coordination number and the composition and symmetry of its immediate environment can be pinpointed.12][13][14][15][16][17][18][19][20][21][22]23 Given this wealth of information, the reasons for presenting another calibration study may not be obvious.Our motivation to calibrate computational Mössbauer spectroscopy is threefold: (i) several technical advances have been made that are not included in previous calibration studies, e.g.newer basis sets, approximations and corrections; (ii) some of the more recent studies do not report on the quadrupole splitting; 19,24 (iii) the emergence of single-atom catalysts, where the catalytically active iron site is embedded in an ill-defined carbonaceous environment with significant π-character, 8 presents a challenging problem for computational Mössbauer spectroscopy and hence warrants a dedicated assessment of its predictive power for this specific coordination environment.
Single-atom catalysts (SACs) are materials synthesized from metal, nitrogen and carbon precursors with at least one pyrolysis step.Typically, the resulting MeNC catalysts are amorphous, carbonaceous materials with multiple phases and contain only small amounts of metal (< 5 wt%) in various forms.3][44] FeNC catalysts show high activity in the oxygen reduction reaction (ORR), a key reaction for fuel cells that are of central importance for green mobility applications.In fact, the ORR activity of recent FeNC catalysts is on par with that of low-platinum content electrodes, the current state of the art. 457][48] An important step towards a better understanding of this promising substitution material for platinum electrodes is to clarify the structure and electronic properties of the active site. 39,43,49,50e current consensus in the literature is that the catalytically active iron ion is surrounded by two to four nitrogen donor atoms embedded in a graphene sheet that may have structural or electronic defects and possibly an additional axial ligand. 39,49,51,52The sketch shown in Figure 1A attempts to summarize the types of environment discussed in the literature.4][55][56][57][58] While other spectroscopies fall short due to the amorphous nature of the SAC material, Mössbauer spectroscopy is ideally equipped to study the coordination environments and electronic structures of iron sites within the amorphous material.However, the definitive assignment of structural characteristics proves difficult purely by comparison with reference data from model complexes.9][60] The long-term goal is to decipher the composition of the active site in its resting state and the changes it undergoes during catalysis in a joint effort of experimental and computational Mössbauer spectroscopy.A Mössbauer spectrum of a typical FeN 4 environment shows a doublet with two defining features: the isomer shift δ and the quadrupole splitting ΔE Q , both sketched in Figure 1B. 1 These arise from the interaction between the charge densities of the iron nucleus and the surrounding electrons.Specifically, the isomer shift reflects the electron density at the iron nucleus, also known as the contact density, and the quadrupole splitting is indicative of an electric field gradient, i.e. the degree of asymmetry within the electron density.While the isomer shift provides information about the iron oxidation state and spin state, the quadrupole splitting can help to differentiate between different electronic states with the same multiplicity.We note that the simultaneous evaluation of both parameters is important to distinguish signals of sites for which the isomer shift or quadrupole splitting by themselves cannot provide an unambiguous assignment (see below).
The numerical values for isomer shift and quadrupole splitting expected for iron in various oxidation and spin states are illustrated in Figure 2.While generally, the isomer shift is higher for lower oxidation states, it can be seen that the observed regions overlap for different oxidation states and different spin states.Fe(II) as one of the relevant oxidation state for the resting state of FeNC-catalysts shows good differentiability between its high spin (S = 2,δ = 0.59-1.45mm s -1 ) and low spin (S = 0, δ = -0.16-0.50mms) states, although intermediate spin (S = 1,δ = 0.26-0.49mm s -1 ) centers are found at the high end of the range expected for low spin complexes.Iron in oxidation state +III is found between -0.17-0.67 mm s -1 with some overlap in the observed regions for all three spin states S = 1/2, S = 3/2 andS = 5/2.Since Fe(II) and Fe(III) may both be present in FeNC catalysts, it is important to note that the isomer shift regions of all ferric spin states overlap with those of ferrous low spin and intermediate spin to some extent.This implies that additional information, such as the quadrupole splitting values shown in Figure 2B, will be needed to assign oxidation and spin states.Lower and higher oxidation states are shown for completeness; while Fe(I) is less relevant to the FeNC intermediates during ORR, higher oxidation states are likely important for the later stages of catalysis.
0][21][22] Details on the computational approach and the expected error margins are given below.A recurring problem in computational iron chemistry is the prediction of the correct spin state energies, 61 the correct spin state being obviously very important to achieve a reliable Mössbauer prediction.2][63] In other cases of course, the electronic structure will have non-negligible multireference character or substantial mixing of low-lying excited states, and in such scenarios DFT is prone to failure.Wavefunction approaches such as the complete active space SCF or density matrix renormalisation group methods for large active spaces in combination with extensive basis sets can lead to accurate predictions of relative spin state energies. 64,65The more recently introduced multiconfiguration pair-density functional theory also shows promising results. 66,67While one might thus think that post-Hartree-Fock-methods are the ideal approach to obtain both more accurate spin state energetics and more reliable Mössbauer parameter predictions, these types of calculation remain far from routine for large molecules with complicated electronic structures.For the problem targeteda screening of iron sites embedded in extended π-systems, paralleling the strategy that has contributed greatly to resolving many questions in bioinorganic chemistry 13,15,20,23,68,69 -a density functional theory based approach is clearly better suited.Given the limitations and challenges set out above, it is however equally clear that the reliability and pitfalls of Mössbauer parameter predictions must be carefully calibrated specifically for FeN 4 sites.In this contribution, computational Mössbauer spectroscopy is calibrated using a set of 20 complexes with FeN 4-6 environments, specifically chosen such that they are representative of plausible active site structures in FeNC catalysts.We find that the TPSSh, B3LYP and PBE0 density functionals are able to predict isomer shift and quadrupole splitting values with approximately equal accuracy.The mean absolute and maximum deviations for the isomer shift are ca.0.05-0.06mm s -1 and 0.12-0.13mm s -1 , and for the quadrupole splitting mean absolute and maximum deviations of ca.0.23-0.26mm s -1 and 0.57-0.93mm s -1 are found.Important raw data including the computed densities is provided in the SI for future use in gauging deviations between computational calibration studies.In addition, we introduce an interactive notebook that based on the reference data set presented here derives predicted Mossbauer parameters with individual associated error margins from computed contact densities and quadrupole splitting values.Defining the computational trust region as twice the mean absolute deviation obtained in the calibration study, we discuss for which of the Mossbauer signatures typically observed in FeNC catalysts one can expect that accurate computational models will be able to differentiate them.

Overview
In Mossbauer spectroscopy, γ-radiation is absorbed by a nucleus. 1 For an absorption event to be detectable, it needs to occur without recoil, which can be significant due to the high energies of the γ-photons (14.4 keV in the case of 57 Fe).Unlike in gaseous or liquid samples, the recoil energy in solid materials can be taken up by the environment of the absorbing atom, i.e. it is transferred into vibrational degrees of freedom of the molecule or crystal.The fraction of recoil-free absorption events is given by the Lamb-Mössbauer factor f where E γ is the energy of the absorbed γ-photon, <x 2 > is the mean square displacement of the nucleus from its equilibrium position, is Planck's constant divided by 2π, and c is the speed of light.Thef -factor is relevant to the total intensity of the recorded absorption.The measurement of Mössbauer spectra relies on the ingenious application of the Doppler effect to achieve partial to full overlap between the emission and absorption lines of the γ-source and the sample, respectively.Details are given in Ref. 1 and references therein.
The typical Mössbauer doublet sketched in Figure 1B arises as a normal absorption line shape at the position of the diagnostic isomer shift that is split due to a magnetic quadrupole interaction of characteristic magnitude.Both parameters, isomer shift and quadrupole splitting (see Figure 1B), arise from the hyperfine interaction of nuclear magnetic dipole moment and electron charge distribution.A third characteristic, the electric quadrupole interaction, is less relevant in this context and only mentioned for completeness. 1

Isomer shift
The isomer shift δ is due to the Coulomb interaction between the absorbing nucleus and the surrounding electrons, expressed as showing the dependence of the isomer shift on the electron density of the absorber nucleus, expressed as the square of the wavefunctionΨ(0) at the nucleus A multiplied with the elemental chargee which is included in the constant α. 71 The term C describes the density at the source nucleus, which is approximated as constant (C = |Ψ (0)| 2 source ).In addition to the elemental charge e , the constant α contains the atomic number Z , the radius of the nucleus R , the nuclear transition energy E 0 , and the electric constantε 0 and the speed of light c as natural constants.ΔR /R describes the relative change of the nuclear radius upon excitation, which is determined experimentally or with the help of quantum chemical calculations of the isomer shift. 71Overall, the isomer shift thus depends chiefly on the electronic charge density at the absorber nucleus.While the electron density is the key quantity in computational Mössbauer spectroscopy, the interpretation of Mössbauer spectra relies on chemical concepts that are mostly formulated in terms of molecular orbitals. 14,72ly electrons in core and valence s-orbitals have a finite probability density to exist at the position of the nucleus; a direct influence of the valence molecular orbitals on the isomer shift is therefore only possible via their s-orbital character.An electron in any orbital with a higher angular momentum (p, d, f, etc.) has zero probability density at the position of the nucleus due to these orbitals' nodal planes.Regardless, the higher angular momentum valence electrons influence the isomer shift indirectly due to shielding of the nuclear charge; for instance, higher iron oxidation states will shield the nuclear charge less effectively, leading to a higher s-density at the nucleus and thus a lower isomer shift (see Figure 2A).Note that for 57 Fe, ΔR /R in Eq. ( 2) is negative, leading to an inverse relationship between s-electron density and isomer shift.
Chemical factors that influence the isomer shift include the oxidation state of iron, iron-ligand bond lengths, covalency and nature of its bonds, electronegativity of the ligands and shielding due to the 3d orbital occupation pattern.All of these factors are important and should be evaluated carefully for a complete interpretation of isomer shifts, however it appears futile to attempt to fully disentangle all individual contributions.A thorough discussion of these factors is presented in Ref. 14.
To compute isomer shifts with DFT, the central quantity is the electron density at the point of the nucleus.From a regression analysis of computed electron densities against experimental isomer shifts, fit parameters a and b are extracted according to: The parameter c can be introduced for convenience.With this approach, correlation lines with R 2 -values up to 0.98 21 and maximum deviations of ca.0.1 mm s -1 have been obtained with density functional theory. 21,22,72he nucleus is approximated as a point charge; studies using finite nuclei models have not shown any significant improvement.The question of how the electron density is obtained 14,72 has been detailed in textbooks and dedicated reviews and therefore we will briefly comment only on two points central to computational Mössbauer spectroscopy: the choice of basis set and the choice of relativistic corrections.
Clearly, the basis set needs to be sufficiently large to adequately describe the contact density, where a cusp in the electron density will occur that is inherently difficult to describe with Gaussian basis sets.From previous calibration studies, 14,[19][20][21][22]24 triple-ζ basis sets such as def2-TZVP have shown good agreement with experiment at reasonable computational cost. A ppular choice for iron is the CP(PPP) basis set 14 that was developed specifically for the description of core properties.A direct comparison of an all-element def2-TZVP description vs. CP(PPP) on iron and def2-TZVP on all other elements shows that standard deviations and R 2 values are overall better with the CP(PPP) basis set, albeit at a somewhat increased computational cost.19,22,72 These basis set choices will not reproduce the absolute electron density at the nuclear position; however although they do yield significant deviations from the experimental value these are highly systematic and result in a satisfactory correlation with experiment.14 Relativistic effects have been shown to vary little for different electronic configurations allowing the introduction of a constant scaling factor (not included in Equation 2).This has been discussed extensively in the literature.14 It was furthermore shown that there is negligible variation of 1s and 2s electron densities as compared with the small variation in electron densities assigned to 3s orbitals and the substantial variation of 4s-like electron densities.14 These data were also used to explain the success of neglecting scalar relativistic effects, which have little influence on the description of the valence orbitals.14 Instead of a calibration of the computed contact density, Filatov 73,74 has presented an approach for the calculation of isomer shifts from first principles, which is mentioned only briefly for completeness.Using a finite nucleus model, the isomer shift is expressed using the derivative of the electronic energy with respect to the radius of the finite nucleus.Importantly, relativistic and electron correlation effects are incorporated immediately.75 This approach was used to determine α in Eq. 2. 76 An important point in the comparison of experimental and computational isomer shifts is its temperature dependence, which emerges from equation (1): higher temperatures result in larger motions of the nucleus, hence larger values for <x 2 > and lower f -factors.1 Note that <x 2 > and thus f can show anisotropy.71 With increasing temperatures, the second-order Doppler shift appears due to significant thermal motions of the source and absorber nuclei which lead to a relativistic shift in the γ-photon proportional to their mean square velocity.1 While at temperatures of up to 77 K, this effect usually contributes less than -0.02 mm s -1 , the influence at room temperature can be on the order of -0.1 mm s -1 , i.e. on par with common deviations between experiment and DFT prediction.Additionally, low-lying excited states with modified electronic structures and hence different Mössbauer parameters may become significantly populated at higher temperatures. Afair comparison between computational and experimental isomer shifts is therefore guaranteed only at temperatures of a few Kelvin.20 Since computational uncertainties will realistically exceed the influence of the temperature, comparison with experimental data obtained below 80 K appears reasonable.

Quadrupole Splitting
The quadrupole splitting ΔE Q occurs when magnetic interactions between nuclear quadrupole moment and electric field gradient at the nuclear position are present: Here, Q is the nuclear electric quadrupole moment for the nuclearI =3/2 state and V i are the eigenvalues of the tensor representing the environmental electric field gradient, 20 which arises when the field at the nuclear position is inhomogeneous due to deviations of the valence electron distribution from cubic symmetry. 72e sign of ΔE Q depends on the relative energy of the magnetically split substates of the nuclear excited state; in the case of Fe, these are states with I z = ±1/2 andI z = ±3/2.Lippard and coworkers have noted that due to the convention that V 3 be the largest eigenvalue, the sign of the quadrupole splitting can be predicted incorrectly in cases where V 1 is small andV 2 and V 3 are very close in magnitude. 20On the other hand, Pápai and Vankó later showed in their extensive correlation study that the sign is predicted correctly in all instances where it is known experimentally (19 of 66 complexes). 21In an MO picture, a perfectly symmetrical t 2g 3 configuration will produce no quadrupole splitting, while any asymmetry in the ligand sphere and, in the case of ionic species, counter ions will lead to an increase in quadrupole splitting assigned as lattice contributions.The nature of the ligand influences the relative magnitude and orientation of V 1-3 , and thereby the sign of ΔE Q .
The variations in the valence electronic structure are quite subtle, and although evidently they can be measured experimentally, DFT calculations often appear not quite sensitive enough to represent the finer nuances of the asymmetry in electron density.The prediction of quadrupole splittings are therefore associated with larger errors than obtained for the isomer shift: correlation lines with R 2 -values of ca.0.95 and mean absolute errors of 0.22 mm s -1 have been obtained with DFT. 21Furthermore, the quadrupole splitting value can be much more sensitive to changes in temperature, rendering low-temperature measurements important for an adequate comparison between experiment and theory.
Choice of electronic structure method From the above discussion, one may wonder whether complete active space methods are better suited to represent the relevant electronic configurations, as they are in principle able to represent multireference cases and low-lying excited states.One study in this direction is briefly mentioned, although as stated in the introduction, active space approaches will have prohibitive computational cost for their applications in screening problems such as that for the active site of FeNC catalysts.Sadoc et al. 77 used CASSCF/CASPT2 wavefunctions to predict Mössbauer isomer shifts.They found a good correlation between the isomer shift and ρ (0) as the sum of natural orbital densities at the point of the nucleus multiplied with the natural orbital occupation numbers (R 2 =0.984).No direct correlation between the effective d-electron count and the isomer shift was found, in line with previous assessments based on density functional theory. 14Instead, the isomer shift correlated with the covalency of the Fe-ligand bonds -taken as the difference between the formal and the computed charge of the iron ion-and was interpreted as a measure of the deviation from the ionic model.

Methods
All calculations were performed with ORCA suite of programs, version 4.2.1. 78The geometries were optimised from the closest available crystal structure using unrestricted Kohn-Sham DFT calculations with the TPSS density functional. 79Ahlrich's def2-TZVP basis set was used for all elements except for carbon and hydrogen, for which def2-SVP was used. 802][83] Grimme's dispersion correction with Becke-Johnson damping (D3BJ) was employed. 84,85Solvation effects were included with the SMD model, 86 choosing water as the modelled solvent.The grid was set to 6 and integration accuracy was increased to 6.0 in ORCA nomenclature.The convergence criteria for the SCF and the geometry optimisation were set to "tight" in ORCA nomenclature; the only exception to this was complex 8 which was optimised with the "NormalOpt" keyword due to its large size.
For the single point calculations, the hybrid functionals TPSSh, 87 PBE0 88 and B3LYP 89 90 were used in conjunction with the RIJCOSX 82 approximation and a GridX choice of 8. Otherwise the same settings as for the geometry optimisations were used, except for a change in basis set: the CP(PPP) 91 basis set was used for Fe and def2-TZVP 80 basis set for all other elements.Additionally, the integration grid for iron was increased to 7.0 in ORCA nomenclature.While all reported calculations use the SMD solvation model for water, we note that omission of the solvation model or its inclusion of with dielectric constants ranging from ε = 4-80 altered the Mössbauer parameters very little.
All calibrations and corresponding visualizations were performed with OriginPro 9.1, 92 except for those related to bootstrapping, which were performed with Python 3.6.8 93(including the packages Numpy 1.18.1, 94 Scikit-learn 0.22.2, 95Matplotlib 3.1.1, 96and Seaborn 0.9.0 97 ).Contact densities, fit parameters, and details on the bootstrapping analysis are reported in the Supporting Information.We also make available a free interactive notebook that allows researchers to request predictions of Mössbauer parameters, including virtual measurement uncertainties (error bars).The notebook can be accessed attinyurl.com/mbs-notebook.

Choice of Reference Systems
We have chosen 20 mononuclear complexes that are representative of the likely iron coordination environment in FeNC catalysts (see Figure 3) and for which Mössbauer data recorded at low temperatures ([?]80 K) are available.This set of complexes contains tetra-, penta-and hexacoordinate iron in oxidation states II, III, IV, V and VI with high spin, intermediate spin and low spin electronic configurations where possible (see Table 1).The ligand spheres are predominantly composed of N-donor ligands, where in many cases the nitrogen atom is part of a π-system such as a porphyrin, bipyridine or phenanthroline ligand.Other ligands are small molecules expected to have σ-and π-interactions with the metal center, e.g.CO, O 2-or N 3-.Ligands that would lead to significant complications in the description of the electronic structure, e.g.nitrosyl ligands, have been excluded.Similarly, no oligonuclear iron complexes have been included to avoid the description of electronic structures that are dominated by antiferromagnetic coupling, e.g. in oxo-bridged dimers, with broken-symmetry DFT.
Although all of the complexes included here except complex 20have been used in previous calibration studies, 21 we note that in the present study none of the ligand spheres were truncated.Hence, the full steric and electronic effects are recovered in the geometry optimizations and subsequent property calculations.

Predicted Mössbauer Parameters
From previous studies it is clear that hybrid functionals are best suited for the prediction of Mössbauer parameters.GGA and meta-GGA functionals perform less favourably, and double-hybrid functionals do not show any significant improvement that would justify the higher computational cost. 19Therefore, three hybrid density functionals, TPSSh (10% HF exchange), B3LYP (20% HF exchange) and PBE0 (25% HF exchange), in combination with the CP(PPP) basis set for iron and the def2-TZVP basis set for all other elements were selected for this calibration study (see Computational Details).We reiterate that the main purpose of this study is to define uncertainty estimates in isomer shift and quadrupole splitting predictions specifically for FeN 4 -environments.S9).Table 3. Experimental isomer shift (δ, mm s -1 ) and quadrupole splitting (ΔE Q , mm s -1 ) values for the 20 complexes along Table 3. Experimental isomer shift (δ, mm s -1 ) and quadrupole splitting (ΔE Q , mm s -1 ) values for the 20 complexes along The predicted values are listed in Table 3; Figure 4 shows that all density functionals perform almost equally well in the prediction of isomer shifts (left) and quadrupole splitting values (right).Note that the calibration lines shown in the main text are obtained after exclusion of complexes 9 and 10 (shown as solid symbols in Figure 4) as is explained in detail below.The fits including all 20 complexes are shown in the Supporting Information (Figure S2).
Beginning with a closer inspection of the isomer shift predictions, the experimental values range from -0.29 mm s -1 to +1.11 mm s -1 .The linear fits obtained with the three density functionals show almost identical R-values of 0.989 (TPSSh), 0.985 (B3LYP) and 0.988 (PBE0, see Table 4).With the B3LYP density functional, the largest discrepancy with experiment is seen for complex9 (0.159 mm s -1 absolute deviation).The maximum absolute discrepancy from the regression line for the B3LYP functional is 0.131 mm s -1 (see Table 4).
For the quadrupole splittings, the experimental values lie between 0.23 mm s -1 and 4.25 mm s -1 in absolute numbers, with the only negative sign reported for complex 17(-1.76mm s -1 ).Since in previous studies the sign was calculated correctly in all cases where it was measured explicitly, 21 the correlation line assumes the sign of the experimental values to match that predicted computationally.As expected, the scatter about the regression line is greater than for the isomer shifts.Overall, the quality of the fit is quite satisfactory with R-values of 0.989 (TPSSh), 0.991 (B3LYP) and 0.987 (PBE0).It can be seen readily that the most significant outliers are complexes9 and 10 , with absolute deviations of 1.601 mm s -1 (9 , B3LYP) and 1.848 mm s -1 (10 , B3LYP) from the experimental value.These much larger than usual deviations led us to exclude them from the fitting procedure, and the underlying electronic structure reasons will be discussed below.Among the set of 18 complexes the fit is based on, the largest deviation is observed for complexes 6 and7 with the B3LYP density functional (0.850 mm s -1 and 0.727 mm s -1 absolute deviation from experiment).The maximum absolute deviation from the B3LYP regression line is smaller at 0.570 mm s -1 .
To summarize, all three density functionals yield roughly the same quality of fit.Because B3LYP provides better predictions for the quadrupole splitting values in terms of the mean absolute deviation and maximum deviations while also performing very well for the isomer shift, any specific discussion in the following makes use of the B3LYP data.We define a trust region for the later evaluation of FeNC catalyst models as twice the mean absolute deviation, i.e. 0.131 mm s -1 for the isomer shift and 0.451 mm s -1 for the quadrupole splitting using the B3LYP density functional.

Bootstrapping Analysis
To examine the reliability of the above calibrations, we conducted a (nonparametric) bootstrapping analysis. 126,127As shown previously with Mössbauer isomer shifts as an exemplary data set, bootstrapping increases the robustness of statistical measures such as fit parameters and relative performances of density functionals. 24Here, we applied Bayesian bootstrapping, 128 which yields smoother results than its original variant. 129The results of the bootstrapping analysis applied to the B3LYP contact density are shown in Figure 5, details are given in the Supporting Information.
The ensembles of regression lines (blue) were obtained by bootstrapping samples from the data set and regressing each sample.The mean over each regression ensemble is marked as a black line and used to make predictions for the isomer shift.The transparent red bands represent 1.96 times the prediction uncertainty (assuming a normal distribution), i.e. it is estimated that 95% of the population is located inside the bands.The results shown in Figure 5 (left) were obtained by bootstrapping all data points except 9 and 10 (not shown) as discussed above.100% of the data lies within the uncertainty band; since each of the 18 remaining data points makes up >5% of the data set, we do not consider this finding a violation of the statistical hypothesis (95% confidence).

Hosted file
image5.emf available at https://authorea.com/users/320493/articles/450023-calibration-ofcomputational-m%C3%B6ssbauer-spectroscopy-to-unravel-active-sites-in-fenc-catalysts-forthe-oxygen-reduction-reactionFor small values of the contact density, the ensemble of regression lines (blue) is skewed toward larger values of the isomer shift.This non-normal distribution of function values is illustrated in the inset (Figure 5A) and can be explained by the presence of the cluster of four data points to the very left of the calibration plot.Neglecting these four data points corresponding to complexes 5 , 6 ,7 and 8 changes the distribution of regression lines and their mean significantly (Figure 5B).Both intercept and slope (which are strongly correlated due to the large values of the contact density) decrease, which also decreases the ability to discriminate between predictions.In the absence of additional data points at intermediate to low contact densities, it is difficult to provide a conclusive answer as to which regression line is more reliable.Nevertheless, we consider the cluster of four data points a valuable addition to the data set for two reasons: (i) the complexes associated with this cluster (5 , 6 , 7 , and 8 ) have different structural motifs, providing an argument against a systematic bias; (ii) coefficients of linear regression models (intercept and slope) tend to be biased toward small absolute values ("regression toward the mean" 130 ).
In an effort to future-proof the calibration presented here and the statistical analysis, we constructed a tool to include more data points, facilitating manifest statistical conclusions beyond the data reported here.To this end, an online database is set up (tinyurl.com/mbs-notebook),which is publicly accessible and open to submissions from other researchers.This database can be used in at least three ways: 1. Obtaining a predicted isomer shift or quadrupole splitting including the associated uncertainty estimates simply by typing in thecomputed contact density or quadrupole splitting.2. Submitting reference data points for additional complexes to obtain more rigorous statistics; the data points will be reviewed by the authors.3. Obtaining complete statistical analyses by submitting new data sets computed with a different computational setup, e.g.different basis sets, solvation models, relativistic corrections, etc.; the data sets will be validated by the authors.
With this database, the authors provide a tool for the prediction and rigorous statistical analysis of computed Mössbauer parameters that will hopefully be of value for all researchers interested in the analysis of electronic structures with 57 Fe Mössbauer spectroscopy.

Complicated Electronic Structures
Finding the correct representation of the electronic structures of complexes 5 and 8 , both Fe(II) high spin (S = 2), presented some challenges, which were however fully resolved.Since an Fe(II) high spin species is a likely intermediate or resting state of FeNC catalysts, we discuss these cases briefly.
For complex 5 , the geometry optimisation resulted in a structure with a Mulliken spin population of 3.74 at the iron ion.For the single point calculation however, the initially obtained orbital occupation pattern can be summarized as (xy) 2 (xz) 1 (yz) 1 (z 2 ) 1 (x 2 -y 2 ) 0 with a Mulliken spin population at the iron of 2.94 and the remaining spin delocalised on the ligand as indicated by significant spin population on the carbon atoms (1.01 in total).With this electronic configuration, erroneous Mössbauer parameters of B3LYP = 0.65 mm s -1 ( exp = 1.03 mm s -1 ) and E Q B3LYP = 1.48 mm s -1 (E Q exp = 4.01 mm s -1 ) were found.Of course, the incorrect Mulliken spin populations allowed for a quick identification of the wrong electronic structure description.To resolve this issue, specific orbitals were rotated to yield a more adequate orbital occupation pattern of (xy) 2 (xz) 1 (yz) 1 (z 2 ) 1 (x 2 -y 2 ) 1 with a Mulliken spin population at the iron of 3.77, which was also energetically preferred.Accordingly, the Mössbauer parameters improved to B3LYP = 1.05 mm s -1 ( exp = 1.03 mm s -1 ) and ), i.e. well within the error margins deduced above.
For complex 8 , a similar problem was encountered already at the stage of TPSS geometry optimisation with a Mulliken spin population of 2.77 instead of the expected value close to 4. With this structure, in the B3LYP single point calculation the d(x 2 -y 2 ) orbital is also found to be unoccupied resulting in a Mulliken spin population of 2.91 and inaccurate predictions for the Mössbauer parameters: B3LYP = 0.47 mm s -1 ( exp = 1.05 mm s -1 ) and E Q B3LYP = 2.28 mm s -1 (E Q exp = 4.25 mm s -1 ).After re-optimisation and orbital rotation, an appropriate representation of the Fe(II) high spin electronic structure is found with a Mulliken spin population of 3.80.With this geometric and electronic structure, the Mössbauer parameters improve significantly to B3LYP = 1.06 mm s -1 ( exp = 1.05 mm s -1 ) and E Q B3LYP = 4.30 mm s -1 (E Q exp = 4.25 mm s -1 ).
In Figure 6, the electron densities and spin densities for 5 in both electronic structure variants are shown.
Although the relevant quantity for Mössbauer spectroscopy is the electron density and not the spin density, 72 it can be readily seen that the electron density is not suitable to discuss any electronic structure changes (Figure 6A/C).Spin density plots are sometimes used in the FeNC literature to characterize and discuss the electronic structures. 131In contrast to the indistinguishable electron densities, the spin density plots do show some discernible differences (Figure 6B/D).In the "initial", incorrect electronic structure description where one iron d-orbital was erroneously unoccupied, spin density (yellow) is seen on the nitrogen donor atoms (Figure 6B).Additionally, there is spin density (red) in the "initial" electronic structure descriptions distributed on the ligand, which can be rationalised as a contribution from a ligand-centered orbital that is occupied but unmatched.However, purely by inspection of the spin densities, it is very difficult to ascertain that an adequate electronic structure is obtained.A spin population analysis appears much more suitable, while an analysis of the MO occupation pattern may be even preferred if an accurate description of the quadrupole splitting is of high importance.The outliers excluded from the fit presented above are complexes9 and 10 (see SI for the fit including all 20 complexes), which have been noted as outliers in previous studies as well. 20,21The isomer shifts predicted with B3LYP deviate by 0.077 mm s -1 (9 ) and 0.159 mm s -1 (10 ) from experiment, the respective quadrupole splitting values by 1.601 mm s -1 (9 ) and 1.848 mm s -1 (10 ).While for9 , is thus predicted within the double mean absolute deviation of 0.131 mm s -1 , the calculated value for10 is outside this trust region.For both complexes, the error in the quadrupole splitting prediction is far outside the trust region of 0.451 mm s -1 and the maximum deviation of 0.57 mm s -1 for the other members of the correlation line.In both cases, the iron(II) intermediate spin (2S +1 = 3) ion is found in the square-planar coordination sphere of a decorated porphyrin.The difficulties in describing such Fe(II)-porphyrin electronic structures have been known for many years, and the body of computational chemistry literature pertaining to the correct prediction of the relative spin states is vast. 65To briefly review the key points of discussion, in D 4h symmetry the two lowestlying triplet states are 3 A 2g with a d-orbital occupation pattern of (xy) 2 (z 2 ) 2 (xz) 1 (yz) 1 (x 2 -y 2 ) 0 and 3 E g with degenerate (xy) 2 (z 2 ) 1 (xz) 2 (yz) 1 (x 2 -y 2 ) 0 or (xy) 2 (z 2 ) 1 (xz) 1 (yz) 2 (x 2 -y 2 ) 0 configurations.Additionally, there is a low-lying 5 A 2g state, and it is difficult to predict the energetic gap between the correct triplet ground state and the higher lying triplet and quintet states.However, for the purpose of estimating the uncertainty of Mössbauer parameter predictions, the key point here is not whether the energetic separation of triplet and quintet states is predicted with quantitative accuracy, but how good the quality of the predicted electron density for the specific triplet states is in comparison with experimental data and in relation to the calibration line.
Given that the use of symmetry can have significant effects, 21,132,133 Pápai and Vankó used symmetry constraints to enforce either triplet state in a DFT calculation, leading to quadrupole splittings of 0.40 mm s -1 (10 , 3 A 2g ) and 3.08 mm s -1 (10 , 3 E g ) with the B3LYP functional.Lowering the symmetry to D 2h led mixing of d(xy) and d(z 2 ) orbitals resulting in a computed quadrupole splitting of 1.25 mm s -1 , quite close to the experimental value of 1.51 mm s -1 .The electronic structures in this manuscript were computed without any symmetry constraints, and virtually no mixing of d(xy) and d(z 2 ) character is found, explaining the pronounced deviation from experiment.
While some of the examples discussed in this section may at first glance look almost trivial, it is important to note that an egregious error in the electronic structure can be quickly spotted with the spin population analysis but may be more difficult to identify based on the spin density contour plots.In most cases, as shown above, such a scenario may be resolved with re-optimising the geometry and/or the electronic structure.An error of this magnitude will most likely result in a drastic difference between the experimental and computed isomer shift.In contrast, it is equally possible that while the spin population is correct, the orbital occupation pattern that leads to this spin population is incorrect.This will influence the isomer shift very little or not at all, as was shown in the discussion of complexes9 and 10 ; however, the prediction of the quadrupole splitting will be much more strongly affected as is expected from the fundamental interactions at play in such a scenario.It is thus clear that treating DFT purely like a black-box method is prone to failure, and we hope to have shown with the above examples that a careful inspection of the electronic structures at least at the level of a spin population analysis, but preferably with an analysis of the MO occupation pattern, is mandatory to achieve reliable results.

Implications for FeNC Systems
From the data presented above, the isomer shift can be predicted with reasonable accuracy for any coordination environment represented in the test set.This explicitly includes the square-planar Fe(II) complexes with difficult electronic structures, implying that for the correct oxidation state, multiplicity and coordination geometry the prediction can be good enough regardless of the exact details of the d-electron configuration or precise nature of the spin state.The trust region for isomer shift predictions is taken as double the mean absolute deviation in the correlation plots of Figure 4, i.e. 0.10 mm s -1 (TPSSh), 0.13 mm s -1 (B3LYP) and 0.12 mm s -1 (PBE0) depending on the choice of density functional.In line with results from previous studies, the quadrupole splitting is predicted with lower accuracy.The trust region, i.e. twice the mean absolute deviation among the penta-and tetracoordinate complexes, is 0.51 mm s -1 (TPSSh), 0.45 mm s -1 (B3LYP) and 0.51 mm s -1 (PBE0) for the three density functionals considered.We note that these trust regions should not be considered as absolute values since individual data points have larger deviations from the correlation line.As detailed above, predictions for the quadrupole splitting of square-planar Fe(II) systems are much more sensitive to the valence electronic structure and thus for any computational FeNC model system care must be taken that all plausible electronic configurations are considered.
In the literature, the Mossbauer data of different species are often presented in graphs where the isomer shift and the quadrupole splitting values of the same signal are plotted on the x-and y-axes, respectively.Figure 7A illustrates this concept for the experimental and computed Mossbauer parameters of the 20 reference complexes.It can be clearly seen that deviations between experiment and calculation, highlighted with connecting lines, are significantly more pronounced in the direction of the quadrupole splitting.
Exemplary room temperature Mossbauer parameters observed in FeNC catalysts are shown in Figure 7B in an analogous graph.The selection includes FeNC catalysts prepared from various preparation approaches: starting from porphyrins, 49,[134][135][136][137][138][139][140][141][142] metal organic frameworks, 38,41,143 and polyaniline in combination with an iron source, 38,[144][145][146][147] formation of active sites by an ammonia treatment 38,137,148 and others. 149,150Note that the labels S i correspond to characteristic doublets known as e.g.D1 or D2 in the original references; however, because some distinct signals carried the same or similar labels in different references we have chosen an independent nomenclature here.Signals labelled with roman numerals are obtained with various treatments, e.g.poisoning by sodium sulphite (IV), 142 presence of excess sulphur during the pyrolysis (II and III) or unknown origin 38 of their specific character.

Hosted file
image7.emf available at https://authorea.com/users/320493/articles/450023-calibration-ofcomputational-m%C3%B6ssbauer-spectroscopy-to-unravel-active-sites-in-fenc-catalysts-forthe-oxygen-reduction-reactionThe challenge posed by these systems for the computational chemist is to identify structural and electronic models for the various FeN x sites that will permit the identification of the environment of the catalytically active site.Among the characteristic doublets, S 1 -S 4 fall within a similar range of isomer shifts between 0.28-0.54mm s -1 .The variation within a group of S i signals from different preparations appears approximately as large as the deviation between groups.In contrast, signal S 5 is well separated at an isomer shift of 0.98-1.01mm s -1 .With the B3LYP trust region for the isomer shift of 0.13 mm s -1 , it is evident that S 1 -S 4 will not be distinguishable by different computational models although it can be expected that one will be able to clearly identify S 5 .The quadrupole splitting values of the S 1 -S 4 and S 5 signals are spread over a wider range of 0.71-2.83mm s -1 .This not only allows a better distinction of the individual signals, but also provides clear evidence that different types of local environment are present for the various FeN 4 sites observed experimentally.Taking the B3LYP trust region for the quadrupole splitting of 0.45 mm s -1 , it appears conceivable that S 1 , S 2 and S 3 will be distinguishable while the exclusive consideration of the quadrupole splitting will not enable to differentiate between S2, S4, S5, IV x and III b .If for instance upon poisoning, a chemical connection between one of the S i doublets and a species labelled with a roman numeral in Figure 7B can be made, this could be used as additional information to better characterize the geometric and electronic structures of typical FeN 4 sites.

Conclusions
We have presented a calibration study of three density functionals, TPSSh, B3LYP and PBE0, for their application in computational Mössbauer spectroscopy.All density functionals were found to perform well for both the isomer shift and the quadrupole splitting prediction; R-values exceeded 0.985 in all cases.We defined trust regions as the double mean absolute deviation of the correlation line and note as exemplary values 0.13 mm s -1 for the isomer shift and 0.45 mm s -1 for the quadrupole splitting obtained with the B3LYP density functional.
Besides the notoriously difficult electronic structures of iron(II) intermediate spin (2S +1 = 3) ions in squareplanar coordination geometry, we discuss other cases where an adequate electronic structure description was not found when treating DFT purely as a black-box method.In order to facilitate comparisons of computational data for iron complexes with complicated electronic structures such as those expected in the coordination environment of FeNC catalysts, it is thus strongly recommended to at least report the spin populations on all relevant atoms.The focus of this study was placed on iron environments similar to those thought to be the active sites in FeNC catalysts, thereby guiding the choice of complexes in the reference set.Despite this intentional constrain, the results are very similar to those of Papai et al ., 21 suggesting that the correlation lines obtained here can be used for other systems as well.The data presented here are made available in an online notebook that allows researchers to obtain predicted Mossbauer parameters and the individual uncertainties from the computed values; the notebook (tinyurl.com/mbs-notebook) is open to submission of additional data points for the computational setups presented here and will accept submissions of entire data sets produced with different methodological choices.In this way, a future-proof and ever-growing calibration of computational Mossbauer spectroscopy is provided that ensures a rigorous and directly comparable statistical analysis of different computational approaches.
In native FeNC catalysts, up to five distinct Mossbauer signatures can be found and several other signals can be produced by different treatments.We discuss the trust regions for isomer shift and quadrupole splitting deduced from the calibration study in the context of the relative positions of these experimental values.From this analysis, it appears probable that appropriate computational models will be able to differentiate between these characteristic spectral features if isomer shift and quadrupole splitting are both considered in the analysis.In other words, through a combination of experimental and computational Mossbauer spectroscopy, one will likely be able to identify the structural and electronic basis for the oxygen reduction reaction in FeNC catalysts.

Figure 1 .
Figure 1.(A) Schematic environment of an FeNC site with the iron ion shown as an orange circle.The extent of the graphene-like environment (black lines) is unknown as indicated by grey dashed lines; note that conjugation is not shown.The nitrogen (blue circle) donation may occur from six-or five-membered rings; the latter will result in local distortions and defects (green lines).Axial ligands may be present, but their number and chemical character is unknown (half blue/half red circles).(B) Generic Mössbauer spectrum

Figure 2 .
Figure 2. Typical values for the isomer shift (mm s -1 ) and the absolute value of the quadrupole splitting (mm s -1 ) for iron in oxidation states I-VI (A) and I-IV (B).The reference compounds for the isomer shift are both molecular complexes and solid state materials, whereas for the quadrupole splitting only complexes with an FeN 4 environment were chosen.The values shown in (A) are adapted from Ref. 70 ; the values in (B) are adapted from Ref.8.

Figure 3 .
Figure 3. Graphical representation of the 20 complexes comprising the reference data set.The elements are colour coded as follows; iron: brown, carbon: grey, nitrogen: blue, oxygen: red, sulfur: yellow, chloride: green; hydrogen atoms are omitted for clarity.

Figure 4 .
Figure 4. Correlation of experimental and calculated isomer shifts (left) and quadrupole splitting parameters (right) for the density functionals TPSSh (A, B), B3LYP (C, D) and PBE0 (E, F).The correlation is based on the unfilled symbols; the filled symbols for complexes 9 and10 were excluded.Numerical values are provided in Table 3, fit parameters in Table 4. Correlation plots including all 20 data points and with the densities obtained by Pápai and Vankó 21 are shown in the SI (Figures S1, S2, TableS9).
Applied magnetic field of 6 T. [b] Applied magnetic field of 7 T.

Figure 5 .
Figure 5. Results of a bootstrapping analysis applied to the B3LYP contact density.Regression lines corresponding to bootstrap samples (1000 in each case) are coloured blue, their means are marked as black lines.The red bands represent 95% confidence intervals of prediction (assuming a normal distribution).(A) Result for all data points (i.e.excluding 9 and 10 , not shown).The histogram in the inset shows the (non-normal) distribution of function values for small values of the contact density.For the histogram, 10000 bootstrap samples were drawn.(B) Result for a data set that excludes data points 5 , 6 , 7 , and 8 (coloured gray) from the data set shown in (A).

Figure 6 .
Figure 6.Isosurface plots of the electron density (left) and spin density (right) for complex 5 with initial (A, B) and final (C, D) electronic structure (for complex 8 , see FigureS4in the Supporting Information).The cutoff radii are set to 0.2 for the electron density and 0.005 for the spin density.

Figure 7 .
Figure 7. Isomer shifts and absolute quadrupole splitting values for the 20 reference complexes (A) and selected signals from FeNC catalysts (B).In (A), calculated values are shown as open circles (?) and experimental values as filled squares (); the error bars are taken as the mean absolute deviation.In (B), the labels differ from the nomenclature used in the original papers (S1 corresponds to D1, S2 and S4 to D2-like species, S3 to D3 in porphyrin-type catalysts, S5 to D3 often obtained after an ammonia treatment step or from MOF-based catalysts; roman numerals indicate doublets that originate from a specific treatment or appear only within distinct synthesis routes).The error bars in (B) correspond to standard deviations for the values of different synthesis routes where applicable; for individual data points, e.g.roman numerals, the error bar is taken as the 95 % confidence interval.

Table 2 .
Key structural parameters of the 20 complexes.Distances are given in Å. Equatorial Fe-N distances are given as a #

Table 2 .
Key structural parameters of the 20 complexes.Distances are given in Å. Equatorial Fe-N distances are given as a

Table 4 .
Fit characteristics and statistical data for the regression analyses of isomer shift and quadrupole splitting with the