Multireference Protonation Energetics of a Dimeric Model of Nitrogenase Iron–Sulfur Clusters

Characterizing the electronic structure of the iron–sulfur clusters in nitrogenase is necessary to understand their role in the nitrogen fixation process. One challenging task is to determine the protonation state of the intermediates in the nitrogen fixing cycle. Here, we use a dimeric iron–sulfur model to study relative energies of protonation at C, S, or Fe. Using a composite method based on coupled cluster and density matrix renormalization group energetics, we converge the relative energies of four protonated configurations with respect to basis set and correlation level. We find that accurate relative energies require large basis sets as well as a proper treatment of multireference and relativistic effects. We have also tested ten density functional approximations for these systems. Most of them give large errors in their relative energies. The best performing functional in this system is B3LYP, which gives mean absolute and maximum deviations of only 10 and 13 kJ/mol with respect to our correlated wave function estimates, respectively, comparable to the uncertainty in our correlated estimates. Our work provides benchmark results for the calibration of new approximate electronic structure methods and density functionals for these problems.


I. INTRODUCTION
2][3][4] Extensive biochemical research has revealed that the catalysis in Mo-nitrogenase takes place in the MoFe-protein, which contains a FeMo-cofactor (FeMoco) cluster, with composition MoFe 7 S 9 C(homocitrate), responsible for the N 2 reduction, and a P-cluster, with composition [Fe 8 S 7 Cys 6 ], that transfers electrons to the FeMoco active site. 5,68][9] Given these structures, ab initio electronic structure computation may be applied to determine the binding sites, reaction intermediates and eventually the catalytic mechanism. 10,114][15][16][17][18][19] This is a formidable task owing to the numerous possible binding positions of the added protons 17 and the complicated electronic structure of the FeMo cluster. 20,21All the above calculations have been performed using density functional theory (DFT). 22However, due to the open-shell and multireference nature of the nitrogenase clusters, the reliability of the obtained DFT results has been called into question and the various functionals predict remarkably different results (over 600 kJ/mol a) Electronic mail: hczhai.ok@gmail.comb) Electronic mail: ulf.ryde@compchem.lu.se c) Electronic mail: gkc1000@gmail.comdifference in the predicted stability of different protonation states for the E 4 state). 11,23,24In spite of the large number of open-shell transition metal centers in these clusters, it has been shown that approximate full configuration interaction (FCI) methods, such as the ab initio density matrix renormalization group (DMRG) algorithm, [25][26][27][28][29][30][31][32][33][34][35][36][37] can tackle the qualitative multireference behavior.][40] These studies focused on active space models of the cluster, which were sufficient for a qualitative understanding of the electronic landscape.However, correctly modeling protonation energetics in FeMoco requires calculations that go beyond qualitative accuracy.In particular, quantitative energetics requires a treatment of electron correlation beyond the strongly correlated active space.
In this work, we study the protonation problem in the simpler case of a dimeric iron-sulfur cluster.We compare four representative protonation sites (C, S, Fe or bridging two Fe ions), and try to estimate the lowest energy site and the energetic ordering.We do this within a composite approach where we separately treat multireference electron correlation using DMRG and dynamic correlation beyond the active space using coupled cluster methods.For comparison, we also include several density functional approaches.Overall, we find that multireference effects and correlation in large basis sets are both crucial to describing the protonation energetics, with correlation effects beyond (perturbative) triples being large.Capturing both effects accurately remains challenging within the composite treatment, but we reach sufficient precision to identify the lowest energy protonation site, as well as to order the relative energies of the protonated structures.In contrast, most of the density functionals that we examine yield qualitative and large quantitative errors for this problem.

II. METHODS
We consider the dimer [(SCH 3 ) 2 FeS(CH 2 )Fe(SCH 3 ) 2 ] 4− as a simple model for studying protonation energetics in an iron-sulfur cluster.For this purpose, we added a proton in the following locations: (a) on the bridging CH 2− 2 group, (b) on the bridging S 2− ion, (c) terminally on one Fe ion, or (d) bridging both Fe ions.These locations are representative of potential protonation sites in FeMoco. 17We denote the four protonated configurations HC, HS, HFe and HFe 2 , respectively.Our goal is to identify the lowest-energy structure and to predict the relative energetics of the four structures.Structures.We optimized structures for the four protonated configurations using a broken-symmetry (BS) open-shell singlet ground-state with two antiferromagnetically coupled high-spin Fe(II) centers (net charge -3).Geometry optimization was carried out at the DFT level with the TPSS functional, 41 the def2-SV(P) basis, 42 and using DFT-D3 43 as a dispersion correction.We show the optimized geometries in Fig. 1 and give the coordinates in Supplementary Material Section I.In the optimized structures, we evaluated 2⟨S z ⟩ at each Fe center; these had opposite signs on the different centers and a magnitude ranging from 2.7 to 3.4 depending on the structure.The structures were first designed and used in this work.We summarize the Fe-ligand bond lengths in the four structures in Table I.We see that the Fe-ligand distances vary substantially depending on the protonation site.In particular, protonation of the bridging S 2− or CH 2− 2 ions (thereby decreasing their charge to -1) increases their bond lengths to Fe by ∼0.2 Å. Adding the proton terminally to Fe1 (formally yielding a hydride and Fe(IV)) also increases the bond lengths between this Fe ion and its other ligands.However, when the added proton bridges the two Fe ions, the bond lengths are not much changed.To check the effect of the level of theory used for geometry optimization, we compared the optimized geometry for each structure using the def2-SV(P) and def2-TZVP basis sets at the DFT level with the TPSS functional plus dispersion correction (without the X2C relativistic correction).The largest geometric difference in the bond lengths listed in Table I between the structures optimized with the two bases is 0.04, 0.09, 0.24, and 0.03 Å for the HC, HS, HFe and HFe 2 structures.The large difference for HFe is because there is a change in the general structure to another local minimum (the hydride ion and one SCH 3 ligand change places; this local minimum can also be found with the smaller basis set and then the maximum difference in bond distances is only 0.04 Å).The changes in the energies relative to HC are 16, 17 and 13 kJ/mol for HS, HFe and HFe 2 , respectively (6 kJ/mol for HFe, when compared to the same local minimum).As these differences are much smaller than the energy difference associated with the choice of protonation site (as discussed later), and in the relevant biological context, the geometries will likely change during the course of dynamics, we will not put too much attention on the exact geometry but simply use the def2-SV(P) optimized geometry in this work.Composite approach.For these four structures, we employed a composite energy approach using coupled cluster with singles and doubles (CCSD) and with perturbative triples [CCSD(T)] to estimate dynamic corre-lation 44,45 and DMRG to estimate multireference correlation.Additional corrections were then included for basis-set completeness and relativistic effects.The final composite energy was computed as where E CCSD(T)-TZ is the CCSD(T) energy obtained with the cc-pVTZ-DK basis set, (E DMRG-act − E CCSD(T)-act ) is a multireference correction (estimated in an active space, which will be defined in the following), and ∆E CCSD(T)-CBS is a basis-set correction (specified below).
The mean-field (HF and DFT) and single-reference post-HF calculations were performed in PySCF. 46,47pin-adapted DMRG calculations were performed in block2. 37For each of these terms, in addition to obtaining energies for the composite energy formula, we carried out additional calculations to understand the impact of different approximations and to estimate the reliability of the corrections.We describe the different terms and these aspects below.Coupled cluster calculations.For the coupled cluster (CC) calculations, we started from BS unrestricted reference determinants.To understand the influence of orbitals and importance of triples, we first carried out calculations using the small cc-pVDZ-DK basis set and the exact two-component (X2C) scalar relativistic Hamiltonian [48][49][50] (larger basis set CC calculations are discussed in the basis-set correction section below) and using 40 frozen core orbitals to reduce computational cost.To examine the impact of the orbital choice, we used both Kohn-Sham DFT (with the TPSS or B3LYP functionals [51][52][53] ), as well as Hartree-Fock Slater determinants.In the unrestricted mean-field calculations, we targeted the projected S z = 0 BS state using an initial guess where the spins in the two Fe atoms were coupled antiferromagnetically. 54The expectation value of the total ⟨S 2 ⟩ in the mean-field state ranged from 3.1 to 5.0 for the structures in this work.Starting from these states, we then computed unrestricted CCSD and CCSD(T) energies. 44,45For the DFT reference determinants, we computed CCSD(T) results based on the semi-canonicalized orbitals [which only diagonalize the occupied-occupied and the virtual-virtual blocks in the Fock matrix.This fixes the definition of the (T) correction].
In addition to the low-spin BS state, we also computed the mean-field and CC energies of the high spin (HS) state with S z = 4.With this, we estimated the energy of the pure-spin (PS) singlet state (S = S z = 0) from the Yamaguchi formula 54,55 , where J is the exchange coupling.For simplicity, we used ⟨S 2 PS ⟩ = 0 and ⟨S 2 BS ⟩ and ⟨S 2 HS ⟩ computed at the CCSD level in the above formula, for computing J from both CCSD and CCSD(T) energies.The difference between the BS and PS state can be taken as an estimate of the missing multireference correlation energy arising from spin-recoupling of the Fe centers, but does not capture other types of multireference correlation.DMRG multireference correction.To better estimate the multireference correction, we constructed an active space for a DMRG calculation.We started from a set of (restricted) natural orbitals obtained by diagonalizing the spin-averaged one-particle density matrix (1PDM) of the CCSD wave function calculated above using the cc-pVDZ-DK basis, and then selected orbitals with the occupancy furthest from 0 or 2 as the active space.
Using this active space, we carried out complete active space configuration interaction (CASCI) spinadapted DMRG 31 calculations, computing the PS singlet state (S = 0) energies.Before performing DMRG, we split-localized the nearly doubly occupied orbitals, nearly empty orbitals, and other orbitals, using the Pipek-Mezey localization algorithm. 56The orbitals were then reordered using the Fiedler algorithm. 32The maximum bond dimension in the DMRG calculations was 5000 [SU(2) multiplets].We used a reverse schedule to generate data for DMRG energy extrapolation, and the DMRG extrapolation error was estimated as one fifth of the energy difference between the extrapolated DMRG energy and the DMRG energy computed at the largest bond dimension 32 (see Supplementary Material Section II).For analysis, we also extracted the largest configuration state function (CSF) coefficient from the DMRG wave function, using a deterministic depth-first search algorithm. 57o obtain a multireference energy correction, we also computed the CCSD and CCSD(T) energies in the same active space, using a BS Hartree-Fock reference.The initial guess for the active space BS UHF density matrix was obtained by projecting the BS UHF density matrix in the full space into the active space.The correction was then computed as ∆E DMRG-act = E DMRG-act − E CCSD(T)-act .
To validate the size of the correction, we considered three different active spaces, to verify that the multireference effects were converged.The dimer model in the cc-pVDZ-DK basis contains 180 electrons in 321 spatial orbitals.We constructed the active spaces from the UHF/CCSD natural orbitals (see Supplementary Material Section III): one with 36 orbitals and 48 electrons (36o, 48e), one with 55 orbitals and 48 electrons (55o, 48e) and one with 63 orbitals and 64 electrons (63o, 64e).The uncertainty in the multireference contribution to the relative energies was then estimated crudely as one half the amount of the DMRG multireference correction for the energy difference between HC and HFe 2 (the least and most multireference structures); here, the one half factor is used to be conservative, as it is the largest estimate of the error still compatible with an assumption that the correction improves the result.We further assume that the DMRG multireference correction computed in the cc-pVDZ-DK basis can be used to correct the CCSD(T) relative energies in the complete basis set (CBS) limit, as shown in Eq. (1).Basis-set correction and relativistic contribution.To estimate the CBS limit, we used energies computed in several bases: the UHF energy in cc-pVDZ/TZ/QZ-DK bases, as well as CCSD and CCSD(T) energies in the cc-pVDZ/TZ-DK bases.To estimate the error in the CBS extrapolation, we additionally computed secondorder Møller-Plesset perturbation theory 58,59 (MP2) energies in the cc-pVDZ/TZ/QZ-DK bases.To independently analyze the size of the relativistic correction we also computed the CCSD(T)/def2-SV(P) energies with and without using the X2C Hamiltonian.The CCSD and CCSD(T) calculations were performed with 40 frozen core orbitals (i.e.excluding the 3s and 3p semicore on Fe).
Using the CC correlation energies in the cc-pVDZ-DK and cc-pVTZ-DK bases, we extrapolated to the CBS limit energies using the two-point formula 60 where X = 2 (DZ), Y = 3 (TZ) and taking β = 2.4.For the corresponding mean-field energy at the CBS limit we simply used E ∞ UHF = E QZ UHF .To estimate the error in the CCSD(T) relative energies in the CBS limit, we performed an independent extrapolation with the MP2 energies, and took half of the difference between the DZ/TZ extrapolation and TZ/QZ extrapolation (i.e. the difference from the average of the two) using the same CBS formula with X = 3 and Y = 4) for the MP2 energies, namely where ∆E is the difference in the CCSD(T) energies between HC and HFe (the structures showing the largest difference in the extrapolated MP2 energies), estimated from the extrapolation based on DZ and TZ bases to the CBS limit.
We briefly note that we did not use multireference dynamic correlation methods (such as DMRG-secondorder N -electron valence perturbation theory (DMRG-NEVPT2)) in this work because the different configurations considered in this work have different bonding topologies.This makes it hard to choose a consistent active space that is also small enough to be used with DMRG-NEVPT2.DFT comparisons.For comparison, we computed BS-DFT energies (without the spin-state corrections) using the X2C Hamiltonian and the TPSS, 41 BLYP, 51,52 PBE, 61 B97-D, 62 r 2 SCAN, 63 TPSSh, 64 B3LYP*, 65 B3LYP, 51-53 PBE0, 66 and M06 67 functionals, with the cc-pVQZ-DK basis set and with dispersion corrections from the DFT-D3 method. 43olvation.To estimate the effect of solvation, we additionally computed single-point DFT energies using the cc-pVDZ-DK basis, the TPSS functional, the DFT-D3 dispersion correction, the X2C correction, and the domain-decomposition COSMO solvation model 68 with a dielectric constant ϵ = 4.0.We compared the relative energies with and without the solvation model.We find that solvation greatly stabilizes the negative charges in the model, reducing the number of formally unbound occupied spatial orbitals (i.e. with positive eigenvalues) from 24-25 to less than three, for the four structures.Solvation also leads to a modest, but non-uniform change in the relative energies of the structures (with respect to HC), from +2.3 kJ/mol for E HS − E HC to −16.8 kJ/mol for E HFe − E HC and −27.1 kJ/mol for E HFe2 − E HC .Clearly solvation is important for accurate studies of the biological system.However, its effects can generally be decoupled from that of the correlation level, and thus for the current benchmark study we will henceforth ignore the effects of solvation for simplicity.

III. RESULTS AND DISCUSSION
In Section III A, we first discuss the CC energies with the cc-pVDZ-DK basis.This will allow us to understand some features of correlation in the system, including the influence of orbital choice and the size of triples correction on the relative protonation energies, setting the stage for understanding the reliability of the composite method.In Section III B we discuss the contribution associated with correcting the BS spin states.In Section III C we discuss the detailed multireference corrections entering the composite energy formula from the DMRG and CC calculations.In Section III D we discuss CC calculations in larger basis sets, the CBS extrapolation entering the composite energy, and the size of relativistic effects.We report the final composite energies, the prediction of the lowest energy protonation site and the relative ordering, and the comparison with DFT calculations in Section III E.

A. CC energies: importance of higher-order correlations
We show the energies of the four protonated structures relative to the HC structure from calculations with HF, CCSD and CCSD(T) methods with the cc-pVDZ-DK basis in Table II and Fig. 2. All methods find that the HC structure is the most stable.However, CCSD and CCSD(T) predict a different relative ordering.In addition, there are large quantitative differences in the relative energies, particularly for the HFe and HFe 2 structures.For example, the relative energy of the HS structure decreases by 69 kJ/mol on moving from UHF to UHF/CCSD(T), but that of HFe 2 de-  creases by 216 kJ/mol.Comparing the energy differences from UHF/CCSD and UHF/CCSD(T), we see a sizable contribution from (T) to the absolute and relative energies.Specifically, the absolute (T) corrections for the HC, HS, HFe and HFe 2 structures are −214, −224, −269 and −287 kJ/mol, meaning that relative to the HC structure, HFe 2 is further stabilized by triples by as much as 73 kJ/mol (see Table II).Consistent with this, the CCSD relative energies are observed to be sensitive to the choice of reference orbitals in HFe and HFe 2 .The large (T) corrections to the relative energies highlight the potentially large contribution of higher-order multireference correlations in the relative protonation energies, especially for the H-Fe bond.

B. Spin-state corrections
As the above calculations used a BS reference, part of the missing higher-order correlation could potentially originate from the energy difference between the BS and PS singlet states.In Table III we report the results from the Yamaguchi energy correction to the BS state and the resulting estimate of the PS relative energies.
The PS state correction to the relative energies is shown in the last two lines in Table III.We see that the PS state correction to the relative energies is modest.It is largest for the HFe/HC difference, where it lowers the relative energy by 13 kJ/mol at the CCSD(T) level.Note that as we explicitly compute multireference contributions from DMRG energies below (which are for PS states), we do not use the PS state energy corrections in the composite energy formula.

C. Multireference effects
To obtain a more complete picture of the multireference effects, we next consider ab initio DMRG energies.In Fig. 3(a) and (b) we plot the CCSD and DMRG natural-orbital occupancies in the (36o, 48e) active space.We see that the most fractionally and singly occupied orbitals are all included in the active space, which suggests the active space (and its larger counterparts) should capture the multireference effects in this system.The occupancy patterns of CCSD and DMRG are qualitatively similar.This shows that while (BS) CCSD and CCSD(T) are not usually considered to be multireference methods, they can be qualitatively correct for (spin-averaged) oneparticle quantities, and thus for most conventional anal-TABLE III.Relative single point energies (in kJ/mol) for the BS, high-spin and (estimated) PS singlet states of the protonated Fe dimers computed using different theories with the cc-pVDZ-DK basis and the scalar relativistic X2C Hamiltonian.The energy of the HC structure is used as the reference.

Theory
Energy yses of bonding.The main problem with the energies obtained from the BS CC methods here is the lack of error cancellation for configurations with varying multireference character (namely, the error in the absolute BS-CCSD(T) energy for the HC and HFe structures can be quite different), rather than the complete failure of the CCSD and CCSD(T) methods.
From the DMRG natural-orbital occupancy plot for the PS state (Fig. 3(b)), we see that there are singly occupied orbitals associated with the Fe centers, but additionally zero, one, two and three orbitals with fractional occupancies between 0.2 and 0.8 (or between 1.2 and 1.8; marked by gray in Fig. 3), respectively, for the HC, HS, HFe and HFe 2 structures.This clearly illustrates the trend of increasing multireference character, beyond spin-recoupling of the Fe centers, in this sequence of four structures.
In Fig. 4 we compare the UHF, CCSD, CCSD(T) and DMRG energy differences for the four protonation states and in Table IV we compare the CCSD(T) and DMRG energy corrections to the CCSD relative energies for the individual structures.We see that within the (36o, 48e), (55o, 48e) and (63o, 64e) active spaces, the (T) contri- bution to the energy difference between HC and HFe 2 is 60%, 66% and 63% of the contribution in the full orbital space, respectively.The largest estimated error for the extrapolated DMRG energy (3 kJ/mol) illustrates that the DMRG energies are almost exact on the current scale of the relative energetics.In all cases, the DMRG and (T) corrections to the CCSD relative energies are in the same direction, as indicated by the dashed lines in Fig. 4.
Based on the data above, we can estimate the errors in the CCSD(T) energies for the HC, HS, HFe and HFe 2 structures to be −3, −7, −22, and −32 kJ/mol (from the 36o active space), −6, −10, −19, and −31 kJ/mol (from the 55o active space), or −9, −14, −20, and −28 kJ/mol (from the 63o active space), respectively.The DMRG correction is thus quite small for the HC and HS structures, but larger for the HFe and HFe 2 structures, reflecting the multireference character of the Fe-H bond.From the DMRG bond dimension M = 5000 wave function in the 63o active space, we obtain a largest CSF weight for the four structures of 0.72, 0.66, 0.51 and 0.36, respectively.This confirms that the error in the CCSD(T) energy increases when the multireference character of the structure increases.
In Fig. 5, we show the trends in the correlation effects beyond CCSD as estimated by DMRG (∆E DMRG = E DMRG − E CCSD ) and (T), for three active space sizes.The curves track each other, justifying the possibility to use the composite energy formula.We use the DMRG results in the 63 orbital active space to correct for the missing multireference effect in the CCSD(T) energies.As discussed in the Methods section, we estimate the uncertainty of this correction for the relative energies as half of the correction for the HFe 2 structure (the largest correction), i.e. ±10 kJ/mol.

D. Basis-set correction and relativistic contribution
In order to study the basis-set effects on the CCSD(T) energies, we computed the UHF, MP2, CCSD and CCSD(T) energies for the four protonated structures using also larger basis sets.The results are listed in Table V and plotted in Fig. 6.We can see that the basisset dependence of the UHF and correlation energies is very different, largely depending on whether the proton is bound to the metal or not.For the HC and HS structures, the UHF relative energies increase (become more positive) as the basis-set size increases, while the CCSD contributions decrease; for the HFe and HFe 2 structures, the trends are opposite.As a result, the basis-set dependence of the mean-field and correlation energies partially cancel, and overall the total CCSD(T) relative energies change non-monotonically with increasing basis-set size.The UHF energies converge at the QZ level and the (T) corrections converge at the TZ level.Therefore, the CCSD(T) relative energy basis-set trend beyond TZ (bottom right panel, Fig. 6) is dominated by the basis-set trend of the CCSD relative energies beyond TZ (top left panel, Fig. 6).
Using the difference between the DZ/TZ-and TZ/QZ-CBS extrapolation energies computed at the MP2 level, we estimate the error at the DZ/TZ-CBS CCSD(T) level to be ±8 kJ/mol for the relative energy of the various structures.TABLE V.The UHF, MP2 correlation, CCSD correlation and (T) correction energies computed using different basis sets.To better highlight trends in the relative energetics, we show the total or correlation energy averaged over the four structures for each basis set.This is then used as a reference energy.It is also interesting to break out the scalar relativistic contributions to the relative energies of the different structures, shown in Table VI.For simplicity we used the def2-SV(P) basis set for a qualitative assessment.We see that the scalar relativistic contribution is important for the relative energies of HFe and HFe 2 [11-12 kJ/mol at the CCSD(T) level].Relativistic effects are clearly necessary to describe differential bonding to Fe.

E. Final composite energies and analysis
In Table VII we summarize our final estimates for the relative energies of the four protonated structures obtained with the composite formula.We show the various contributions to the energy differences in Fig. 7. Overall, we find that E HC < E HS < E HFe2 < E HFe .
Both the basis-set and high-order correlation effects are important to obtain the correct qualitative ordering.While CCSD(T)/TZ may often be considered to produce reasonable results for the thermochemistry of small molecules, this is not the case for the Fe-S clusters: multireference effects beyond (T) and basis-set effects beyond TZ change the relative protonation energy of HC and HFe 2 by −19 and +21 kJ/mol, respectively.As we needed to perform extrapolations to obtain both the multireference and basis-set corrections, our estimated uncertainty in these energies is ±10 and ±8 kJ/mol, respectively.However, it must be stressed that our estimates of the uncertainties are quite crude.Interestingly, although the multireference and basis-set contributions are individually large, they have opposite signs.Consequently, the combined contribution is significantly smaller, and more closely resembles the raw CCSD(T)/TZ result.
Table VII and Fig. 7 also include relative energies calculated with ten different DFT methods.It can be seen that the BLYP, B97-D, r2SCAN, TPSSh, B3LYP*, B3LYP, and PBE0 functionals all obtain the correct qualitative ordering of the structures, while the TPSS, PBE, and M06 functionals do not.Out of the functionals with the correct ordering, the standard hybrid functionals B3LYP and PBE0 recover the composite method energetics to (approximately) within the estimated uncertainty in our composite results (18 kJ/mol, from adding the uncertainty in the multireference and basis-set correction) while the other functionals do not.Overall, there is a wide spread in the DFT predictions, for example, the range of the energy difference between HC and HS differs from our best estimate by 0-45 kJ/mol.The largest errors are found for the HFe 2 structure, where the PBE functional gives an error in the relative protonation energy of 112 kJ/mol.These effects are expected to be multiplied when there are multiple protons involved, as is the case for the E 4 intermediate state of FeMoco.Our results are thus consistent with the large variation in protonation energies (hundreds of kJ/mol) observed when using different functionals to study multiply-protonated structures in the E 4 intermediate 23 .In this work, we studied the protonation energetics of a dimeric iron-sulfur cluster, as a simple model for the protonation of nitrogenase iron-sulfur clusters, such as the intermediates in FeMoco.Using a composite method based on CC and DMRG energies, we estimated the relative protonation energies of four representative structures (protonated on C, S, Fe or bridging two Fe) in the multireference and basis-set limits.We found that both multireference and basis-set effects are extremely important to capturing the correct energy ordering.Importantly, even though we are studying the seemingly simple process of adding a single proton to the cluster, basis-set effects beyond triple-zeta, and correlation effects beyond (perturbative) triples, contribute about 20 kJ/mol to the relative energies (although the contributions have opposite signs).This highlights the challenge of computing accurate energetics for even larger clusters.
The current work relies on a number of extrapolations to obtain the basis-set and correlation-effect limits.These extrapolations, as well as the crude estimates of the errors associated with them, are not entirely sat- isfactory.While some of these steps could be removed by performing more demanding computations, it may be challenging to scale such a strategy to the larger ironsulfur clusters.In particular, while perturbative triples formed a reasonable starting point for the relative energetics in this cluster, it is unclear whether this will be the case in larger iron-sulfur clusters.The density functionals that we examined yielded a wide range of predictions, from qualitatively incorrect results, to results compatible with our estimates, with the hybrid B3LYP functional giving the best results.The different behaviors of the functionals highlights the well-known importance of tailoring the functional in challenging transition metal problems.Benchmark energetics, such as those from this work, thus serve as a starting point for choosing appropriate functionals to explore the chemistry of larger Fe-S clusters.

SUPPLEMENTARY MATERIAL
Cluster geometries, DMRG energy extrapolations performed in the three active spaces, and the figures of the CCSD natural orbitals for constructing the active space.

FIG. 1 .
FIG. 1.The geometries of the protonated dimer complex [HFe2S(CH2)(SCH3)4] 3− with the added proton on (a) C, (b) S, (c) terminally on one Fe atom and (d) bridging both Fe atoms.The added proton is shown in purple.

FIG. 2 .
FIG. 2. Relative single point energies of the BS state of the protonated Fe dimers computed using different theories with the cc-pVDZ-DK basis and the scalar relativistic X2C Hamiltonian.The energy of the HC struture is used as the reference.

FIG. 3 .
FIG. 3. Natural orbital occupancies computed using (a) CCSD, and (b) DMRG in the (36o, 48e) active space, for the four protonated structures.The gray area shows the range of fractional occupancy, as defined in the main text.

FIG. 4 .
FIG. 4. Relative single point energies for the protonated Fe dimers computed using UHF, CCSD, CCSD(T) and DMRG with the cc-pVDZ-DK basis, in the full orbital space and in the (36o, 48e) and (55o, 48e) active spaces.The energy of the HC structure is used as the reference.The trends in the (T) and DMRG correction for the relative energies are shown by dashed lines.

FIG. 5 .
FIG. 5. CCSD(T) and DMRG correlation energies, relative to CCSD, of the protonated Fe dimers computed for active spaces of different sizes (using UHF/CCSD natural orbitals) in the cc-pVDZ-DK basis set.

5 TZ 3 FIG. 6 .
FIG. 6. Trends in the energies of the protonated Fe dimers for (a) the CCSD correlation energies, (b) the (T) corrections, (c) the UHF energies and (d) the total CCSD(T) energies, as a function of basis set.For each basis set, the total or correlation energies are shifted by their average among the four structures.For UHF and CCSD(T) energies of the HC structure, an additional +100 kJ/mol shift is added for clarity.

FIG. 7 .
FIG. 7. Comparison between the difference in single-point energy of the protonated Fe dimers computed using the composite CCSD(T)/DMRG approach and DFT with different functionals.The energy of the HC structure is used as the reference.The cc-pVQZ-DK basis set and CBS extrapolation results are used for mean-field and post-mean-field methods, respectively, unless otherwise specified.The uncertainty in the energy difference is shown as the error bar.

TABLE I .
The bond lengths (in Å) between the iron ions (Fe1 and Fe2) and their direct ligands in the optimized geometries of the four dimer models.St1-St4 are the sulfur atoms of four terminal SCH − 3 groups, S br is the bridging S 2− ion, C br is the carbon atom of the bridging CH 2−

TABLE II .
Relative single-point energies (in kJ/mol) for the BS state of the four protonated Fe dimer structures computed using different theories with the cc-pVDZ-DK basis and the scalar relativistic X2C Hamiltonian.The energy of the HC structure is used as the reference.

TABLE VI .
Scalar relativistic corrections (in kJ/mol) to the relative energies computed using different theories with the def2-SV(P) basis.The energy of the HC structure is used as the reference energy for all energies.∆E ref represents the energy difference with no relativistic corrections.

TABLE VII .
Relative single-point energies (in kJ/mol) for the four protonated Fe dimer structures computed by the composite method.The energy of the HC structure is used as the reference.UKS energies with different DFT functionals are also listed for comparison.The last two columns show the mean and max deviation in the energy difference E − EHC between the composite and UKS methods for each functional.