# Multireference protonation energetics of a dimeric model of nitrogenase
iron–sulfur clusters

Huanchen Zhai<EMAIL_ADDRESS>Division of Chemistry and Chemical
Engineering, California Institute of Technology, Pasadena, CA 91125, USA
Seunghoon Lee  Division of Chemistry and Chemical Engineering, California
Institute of Technology, Pasadena, CA 91125, USA  Zhi-Hao Cui  Division of
Chemistry and Chemical Engineering, California Institute of Technology,
Pasadena, CA 91125, USA  Lili Cao  Department of Theoretical Chemistry, Lund
University, P. O. Box 124, SE-221 00 Lund, Sweden  Ulf Ryde
<EMAIL_ADDRESS>Department of Theoretical Chemistry, Lund University,
P. O. Box 124, SE-221 00 Lund, Sweden  Garnet Kin-Lic Chan<EMAIL_ADDRESS>Division of Chemistry and Chemical Engineering, California Institute of
Technology, Pasadena, CA 91125, USA

###### Abstract

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 the 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 wavefunction 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

Nitrogenase is the only enzyme that can catalyze the conversion of atmospheric
dinitrogen ($\mathrm{N_{2}}$) to ammonia ($\mathrm{NH_{3}}$), the key reaction
in nitrogen fixation.[1, 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
$\mathrm{MoFe_{7}S_{9}C(homocitrate)}$, responsible for the $\mathrm{N_{2}}$
reduction, and a P-cluster, with composition $\mathrm{[Fe_{8}S_{7}Cys_{6}]}$,
that transfers electrons to the FeMoco active site.[5, 6] During the last two
decades, the atomic structures of the nitrogenase clusters have been
determined by X-ray crystallography.[7, 8, 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, 11]

Recently, the intriguing $\mathrm{E_{4}}$ intermediate of nitrogenase,[12]
formed by adding four electrons and protons to the $\mathrm{E_{0}}$ resting
state of FeMoco and responsible for the binding of N2, has been studied
computationally by several groups.[13, 14, 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,
21] All the above calculations have been performed using density functional
theory (DFT).[22] However, 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 difference in the predicted stability of different
protonation states for the E4 state).[23, 24, 11] In 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. For example, earlier calculations using DMRG provided new insights
into the electronic structure of the P-cluster with its manifold of low-energy
electronic states and their non-classical spin correlations.[38, 39, 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
$[\mathrm{(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 S2- ion, (c) terminally on one Fe
ion, or (d) bridging both Fe ions. These locations are representative of
potential protonation sites in FeMoco.[17] We denote the four protonated
configurations HC, HS, HFe and $\mathrm{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\langle S_{z}\rangle$ 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.

Figure 1: The geometries of the protonated dimer complex
$\mathrm{[HFe_{2}S(CH_{2})(SCH_{3})_{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.

We summarize the Fe–ligand bond lengths in the four structures in Table 1. We
see that the Fe–ligand distances vary substantially depending on the
protonation site. In particular, protonation of the bridging S2- or
CH${}_{2}^{2-}$ ions (thereby decreasing their charge to –1) increases their
bond lengths to Fe by $\sim$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.

Table 1: 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, Sbr is the bridging S2- ion, Cbr is the carbon atom of the bridging CH${}_{2}^{2-}$ ion and H is the added proton, when binding to Fe.  |  |  |  |
---|---|---|---|---
Bond |  HC |  HS |  HFe |  $\mathrm{HFe_{2}}$
|  |  |  |
|  |  |  |
St1$-$ Fe1 |  2.45 |  2.43 |  2.62 |  2.43
St2$-$ Fe1 |  2.45 |  2.45 |  2.50 |  2.48
St3$-$ Fe2 |  2.45 |  2.44 |  2.42 |  2.39
St4$-$ Fe2 |  2.42 |  2.39 |  2.38 |  2.43
|  |  |  |
Sbr$-$ Fe1 |  2.26 |  2.47 |  2.42 |  2.31
Sbr$-$ Fe2 |  2.28 |  2.53 |  2.21 |  2.34
|  |  |  |
Average (all S$-$Fe) |  2.38 |  2.45 |  2.42 |  2.40
|  |  |  |
Cbr$-$Fe1 |  2.20 |  1.98 |  2.02 |  1.97
Cbr$-$Fe2 |  2.19 |  1.99 |  1.95 |  1.99
|  |  |  |
H$-$Fe1 |  |  |  1.76 |  1.72
H$-$Fe2 |  |  |  |  1.76

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 1 between the structures
optimized with the two bases is 0.04, 0.09, 0.24, and 0.03 Å for the HC, HS,
HFe and $\mathrm{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 $\mathrm{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 $\mathrm{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 correlation [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

$E_{\textrm{composite}}=E_{\textrm{CCSD(T)-TZ}}\
+(E_{\textrm{DMRG\mbox{-}act}}-E_{\textrm{CCSD(T)\mbox{-}act}})\\\ +\Delta
E_{\mathrm{CCSD(T)\mbox{-}CBS}}$ ((1))

where $E_{\textrm{CCSD(T)-TZ}}$ is the CCSD(T) energy obtained with the cc-
pVTZ-DK basis set,
$(E_{\textrm{DMRG\mbox{-}act}}-E_{\textrm{CCSD(T)\mbox{-}act}})$ is a
multireference correction (estimated in an active space, which will be defined
in the following), and $\Delta E_{\mathrm{CCSD(T)\mbox{-}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, 47] Spin-adapted DMRG calculations were performed in
block2.[37] For 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.[54] The expectation value of the total
$\langle S^{2}\rangle$ 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, 45] For 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[55, 54]

$J=\frac{E_{\mathrm{HS}}-E_{\mathrm{BS}}}{\langle
S^{2}_{\mathrm{BS}}\rangle-\langle
S^{2}_{\mathrm{HS}}\rangle}=\frac{E_{\mathrm{PS}}-E_{\mathrm{BS}}}{\langle
S^{2}_{\mathrm{BS}}\rangle-\langle S^{2}_{\mathrm{PS}}\rangle},$

where $J$ is the exchange coupling. For simplicity, we used $\langle
S^{2}_{\mathrm{PS}}\rangle=0$ and $\langle S^{2}_{\mathrm{BS}}\rangle$ and
$\langle S^{2}_{\mathrm{HS}}\rangle$ 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) spin-adapted 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.[56] The orbitals were then
reordered using the Fiedler algorithm.[32] The 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.[57]

To 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 $\Delta
E_{\mathrm{DMRG\mbox{-}act}}=E_{\mathrm{DMRG\mbox{-}act}}-E_{\mathrm{CCSD(T)\mbox{-}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 $\mathrm{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 second-
order 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]

$\displaystyle
E_{\mathrm{corr}}^{\infty}=\frac{X^{\beta}E_{\mathrm{corr}}^{(X)}-Y^{\beta}E_{\mathrm{corr}}^{(Y)}}{X^{\beta}-Y^{\beta}}$

where $X=2\,(\mathrm{DZ}),Y=3\,(\mathrm{TZ})$ and taking $\beta=2.4$. For the
corresponding mean-field energy at the CBS limit we simply used
$E_{\mathrm{UHF}}^{\infty}=E_{\mathrm{UHF}}^{\mathrm{QZ}}$.

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

$\Big{|}\Delta E_{\mathrm{CCSD(T)}}^{\mathrm{DZ/TZ}\rightarrow\infty}-\Delta
E_{\mathrm{CCSD(T)}}^{\mathrm{exact}}\Big{|}\\\
\approx\frac{1}{2}\Big{(}\Delta
E_{\mathrm{MP2}}^{\mathrm{TZ/QZ}\rightarrow\infty}-\Delta
E_{\mathrm{MP2}}^{\mathrm{DZ/TZ}\rightarrow\infty}\Big{)}$ ((2))

where $\Delta E_{\mathrm{CCSD(T)}}^{\mathrm{DZ/TZ}\rightarrow\infty}$ 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-second-order $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] r2SCAN,[63] TPSSh,[64] B3LYP*,[65] B3LYP,[51, 52, 53]
PBE0,[66] and M06[67] functionals, with the cc-pVQZ-DK basis set and with
dispersion corrections from the DFT-D3 method.[43]

Solvation. 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 $\epsilon=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\ \mathrm{kJ/mol}$ for
$E_{\mathrm{HS}}-E_{\mathrm{HC}}$ to $-16.8\ \mathrm{kJ/mol}$ for
$E_{\mathrm{HFe}}-E_{\mathrm{HC}}$ and $-27.1\ \mathrm{kJ/mol}$ for
$E_{\mathrm{HFe_{2}}}-E_{\mathrm{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.1, 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.2 we discuss the
contribution associated with correcting the BS spin states. In Section III.3
we discuss the detailed multireference corrections entering the composite
energy formula from the DMRG and CC calculations. In Section III.4 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.5.

###  III.1 CC energies: importance of higher-order correlations

Table 2: 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. Theory | Energy difference $E-E_{\mathrm{HC}}$
---|---
|  HC |  HS |  HFe |  $\mathrm{HFe_{2}}$
|  |  |  |
UHF |  0.0 |  203.8 |  322.5 |  396.3
|  |  |  |
UHF/CCSD |  0.0 |  143.7 |  245.3 |  253.1
UKS-TPSS/CCSD |  0.0 |  154.9 |  278.9 |  305.9
UKS-B3LYP/CCSD |  0.0 |  149.3 |  269.6 |  286.8
|  |  |  |
UHF/CCSD(T) |  0.0 |  134.3 |  190.1 |  179.9
UKS-TPSS/CCSD(T) |  0.0 |  132.6 |  179.4 |  163.8
UKS-B3LYP/CCSD(T) |  0.0 |  132.4 |  185.4 |  167.3
|  |  |  |
Figure 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.

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 2 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 $\mathrm{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 $\mathrm{HFe_{2}}$ decreases 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
$\mathrm{HFe_{2}}$ structures are $-214,-224,-269$ and $-287\
\mathrm{kJ/mol}$, meaning that relative to the HC structure,
$\mathrm{HFe_{2}}$ is further stabilized by triples by as much as 73 kJ/mol
(see Table 2). Consistent with this, the CCSD relative energies are observed
to be sensitive to the choice of reference orbitals in HFe and
$\mathrm{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.

###  III.2 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 3 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 3. 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.

Table 3: 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 difference $E-E_{\mathrm{HC}}$
---|---
|  HC |  HS |  HFe |  $\mathrm{HFe_{2}}$
|  |  |  |
High-spin $S_{z}=4$
|  |  |  |
UHF/CCSD |  0.0 |  164.1 |  291.6 |  309.0
UHF/CCSD(T) |  0.0 |  158.5 |  230.7 |  211.9
CCSD $\langle S^{2}\rangle$ |  20.01 |  20.01 |  20.32 |  20.08
|  |  |  |
BS singlet
|  |  |  |
UHF/CCSD |  0.0 |  143.7 |  245.3 |  253.1
UHF/CCSD(T) |  0.0 |  134.3 |  190.1 |  179.9
CCSD $\langle S^{2}\rangle$ |  3.89 |  3.78 |  4.57 |  4.15
|  |  |  |
Exchange coupling $J$ (estimated, $\mathrm{cm}^{-1}$)
|  |  |  |
UHF/CCSD |  $-94.0$ |  $-198.2$ |  $-342.0$ |  $-388.5$
UHF/CCSD(T) |  $-112.2$ |  $-236.2$ |  $-329.9$ |  $-281.5$
|  |  |  |
PS singlet (estimated)
|  |  |  |
UHF/CCSD |  0.0 |  139.1 |  231.0 |  238.2
UHF/CCSD(T) |  0.0 |  128.8 |  177.3 |  171.1
|  |  |  |
PS singlet correction (estimated)
|  |  |  |
UHF/CCSD |  $-4.4$ |  $-9.0$ |  $-18.7$ |  $-19.3$
UHF/CCSD(T) |  $-5.2$ |  $-10.7$ |  $-18.0$ |  $-14.0$
|  |  |  |

###  III.3 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) one-particle quantities,
and thus for most conventional analyses 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 $\mathrm{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.

Figure 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. Figure 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.  Table 4: Comparison between UHF/CCSD(T) and DMRG energy corrections (in kJ/mol) for individual protonated structures computed using the cc-pVDZ-DK basis. The CCSD energy is used as the reference for CCSD(T) and DMRG energies. Note that the DMRG energies are computed for the PS singlet state while other energies are computed for the BS state. Energies and errors are listed for individual structures unless otherwise specified by “(relative to HC)”. Theory | Energy correction $E-E_{\mathrm{CCSD}}$
---|---
|  HC |  HS |  HFe |  $\mathrm{HFe_{2}}$
|  |  |  |
active space (36o, 48e)
|  |  |  |
CCSD (relative to HC) |  0.0 |  169.2 |  241.8 |  290.7
CCSD(T) |  $-3.4$ |  $-4.4$ |  $-30.4$ |  $-47.0$
DMRG (extrapolated) |  $-6.8$ |  $-11.9$ |  $-52.8$ |  $-79.4$
DMRG extrap. error |  0.0 |  0.0 |  0.1 |  0.2
|  |  |  |
active space (55o, 48e)
|  |  |  |
CCSD (relative to HC) |  0.0 |  160.4 |  232.1 |  269.5
CCSD(T) |  $-10.8$ |  $-13.2$ |  $-42.3$ |  $-59.2$
DMRG (extrapolated) |  $-16.5$ |  $-23.5$ |  $-61.0$ |  $-90.1$
DMRG extrap. error |  0.4 |  0.4 |  1.0 |  3.2
|  |  |  |
active space (63o, 64e)
|  |  |  |
CCSD (relative to HC) |  0.0 |  157.9 |  270.1 |  304.2
CCSD(T) |  $-18.3$ |  $-21.4$ |  $-46.9$ |  $-64.6$
DMRG (extrapolated) |  $-27.1$ |  $-35.5$ |  $-67.2$ |  $-92.6$
DMRG extrap. error |  0.5 |  0.5 |  1.2 |  2.8
|  |  |  |
full orbital space
|  |  |  |
CCSD (relative to HC) |  0.0 |  143.7 |  245.3 |  253.1
CCSD(T) |  $-214.0$ |  $-223.5$ |  $-269.2$ |  $-287.3$
|  |  |  |

In Fig. 4 we compare the UHF, CCSD, CCSD(T) and DMRG energy differences for
the four protonation states and in Table 4 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) contribution to the energy difference between HC and
$\mathrm{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\ \mathrm{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 $\mathrm{HFe_{2}}$ structures to be $-3,-7,-22,$ and
$-32\ \mathrm{kJ/mol}$ (from the 36o active space), $-6,-10,-19,$ and $-31\
\mathrm{kJ/mol}$ (from the 55o active space), or $-9,-14,-20,$ and $-28\
\mathrm{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 $\mathrm{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 ($\Delta
E_{\mathrm{DMRG}}=E_{\mathrm{DMRG}}-E_{\mathrm{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
HFe2 structure (the largest correction), i.e. $\pm 10\ \mathrm{kJ/mol}$.

Figure 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.

###  III.4 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 5 and plotted in
Fig. 6. We can see that the basis-set 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 $\mathrm{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 $\pm 8\ \mathrm{kJ/mol}$ for the relative energy of the various
structures.

Table 5: 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.  Basis |  $E_{\mathrm{average}}$ | $E-E_{\mathrm{average}}$ (kJ/mol)
---|---|---
|  (Hartree) |  HC |  HS |  HFe |  $\mathrm{HFe_{2}}$
|  |  |  |  |
UHF
|  |  |  |  |
def2-SV(P) |  $-4728.8928$ |  $-248.7$ |  $-40.7$ |  $+108.3$ |  $+181.1$
cc-pVDZ-DK |  $-4733.8071$ |  $-230.7$ |  $-26.9$ |  $+91.9$ |  $+165.7$
cc-pVTZ-DK |  $-4733.9598$ |  $-221.5$ |  $-17.9$ |  $+83.2$ |  $+156.1$
cc-pVQZ-DK |  $-4734.0044$ |  $-220.9$ |  $-17.9$ |  $+83.2$ |  $+155.6$
|  |  |  |  |
UHF/MP2 (correlation energy)
|  |  |  |  |
def2-SV(P) |  $-1.6831$ |  $+49.5$ |  $+21.1$ |  $-12.7$ |  $-57.9$
cc-pVDZ-DK |  $-2.3403$ |  $+27.9$ |  $-28.5$ |  $+36.4$ |  $-35.9$
cc-pVTZ-DK |  $-3.3457$ |  $-0.6$ |  $-52.7$ |  $+66.2$ |  $-13.0$
cc-pVQZ-DK |  $-3.9214$ |  $-5.0$ |  $-58.0$ |  $+71.6$ |  $-8.5$
TZ/QZ CBS |  $-4.5001$ |  $-9.5$ |  $-63.3$ |  $+77.0$ |  $-4.1$
DZ/TZ CBS |  $-3.9565$ |  $-17.9$ |  $-67.4$ |  $+84.3$ |  $+0.9$
|  |  |  |  |
UHF/CCSD (correlation energy)
|  |  |  |  |
def2-SV(P) |  $-1.8521$ |  $+81.3$ |  $+39.4$ |  $-34.7$ |  $-86.0$
cc-pVDZ-DK |  $-2.4104$ |  $+70.1$ |  $+10.0$ |  $-7.1$ |  $-73.1$
cc-pVTZ-DK |  $-3.1171$ |  $+48.8$ |  $-5.5$ |  $+13.4$ |  $-56.7$
DZ/TZ CBS |  $-3.5464$ |  $+35.8$ |  $-14.9$ |  $+25.8$ |  $-46.7$
|  |  |  |  |
UHF/CCSD(T) [(T) only]
|  |  |  |  |
def2-SV(P) |  $-0.0701$ |  $+26.8$ |  $+21.1$ |  $-16.9$ |  $-31.1$
cc-pVDZ-DK |  $-0.0946$ |  $+34.5$ |  $+25.0$ |  $-20.7$ |  $-38.8$
cc-pVTZ-DK |  $-0.1491$ |  $+34.6$ |  $+23.8$ |  $-18.7$ |  $-39.7$
DZ/TZ CBS |  $-0.1821$ |  $+34.7$ |  $+23.0$ |  $-17.4$ |  $-40.3$
|  |  |  |  |
Figure 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\ \mathrm{kJ/mol}$ shift is added for clarity.

It is also interesting to break out the scalar relativistic contributions to
the relative energies of the different structures, shown in Table 6. 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 $\mathrm{HFe_{2}}$ [11–12 kJ/mol at the CCSD(T) level].
Relativistic effects are clearly necessary to describe differential bonding to
Fe.

Table 6: 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. $\Delta E_{\mathrm{ref}}$ represents the energy difference with no relativistic corrections. Theory | Energy difference $E-E_{\mathrm{HC}}$
---|---
|  HC |  HS |  HFe |  $\mathrm{HFe_{2}}$
|  |  |  |
Relativistic: $\Delta E_{\mathrm{X2C}}-\Delta E_{\mathrm{ref}}$
|  |  |  |
UHF |  0.0 |  $-3.7$ |  $-21.5$ |  $-21.9$
UKS-TPSS |  0.0 |  $-2.2$ |  $-8.7$ |  $-7.7$
UKS-B3LYP |  0.0 |  $-2.1$ |  $-10.1$ |  $-8.8$
|  |  |  |
UHF/CCSD |  0.0 |  $-2.5$ |  $-13.7$ |  $-13.5$
UHF/CCSD(T) |  0.0 |  $-2.3$ |  $-11.9$ |  $-10.7$
|  |  |  |

###  III.5 Final composite energies and analysis

In Table 7 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_{\mathrm{HC}}<E_{\mathrm{HS}}<E_{\mathrm{HFe_{2}}}<E_{\mathrm{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
$\mathrm{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 $\pm 10$ and $\pm
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 7 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 $\mathrm{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 $\mathrm{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 $\mathrm{E_{4}}$
intermediate [23].

Table 7: 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-E_{\mathrm{HC}}$ between the composite and UKS methods for each functional. Correction/functional | Energy difference $E-E_{\mathrm{HC}}$ | UKS deviation
---|---|---
|  HC |  HS |  HFe |  $\mathrm{HFe_{2}}$ |  mean |  max
|  |  |  |  |  |
UHF/CCSD(T) (with X2C)
|  |  |  |  |  |
uncorrected (cc-pVTZ-DK) |  0.0 |  138.6 |  216.1 |  197.9 |  |
|  |  |  |  |  |
multireference correction |  0.0 |  $-5.2$ |  $-11.4$ |  $-19.2$ |  |
basis-set correction |  0.0 |  $+2.1$ |  $+25.9$ |  $+21.2$ |  |
|  |  |  |  |  |
total |  0.0 |  135.5 |  230.6 |  199.9 |  0.0 |  0.0
|  |  |  |  |  |
UKS (with X2C, DFT-D3 and cc-pVQZ-DK basis)
|  |  |  |  |  |
PBE |  0.0 |  93.1 |  145.6 |  87.8 |  79.8 |  112.1
BLYP |  0.0 |  90.4 |  143.3 |  97.0 |  78.4 |  102.9
TPSS |  0.0 |  98.6 |  161.8 |  97.4 |  69.4 |  102.5
B97-D |  0.0 |  100.2 |  143.0 |  115.2 |  69.2 |  87.6
r2SCAN |  0.0 |  114.5 |  163.5 |  128.3 |  53.2 |  71.6
TPSSh |  0.0 |  115.6 |  200.9 |  156.0 |  31.2 |  43.9
B3LYP* |  0.0 |  113.9 |  203.4 |  176.0 |  24.2 |  27.2
M06 |  0.0 |  135.5 |  190.6 |  194.0 |  15.3 |  40.0
PBE0 |  0.0 |  132.0 |  239.5 |  225.0 |  12.5 |  25.1
B3LYP |  0.0 |  122.1 |  223.8 |  209.1 |  9.8 |  13.4
|  |  |  |  |  |
Figure 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.

##  IV Conclusions

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 satisfactory.
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
iron–sulfur 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.

###### Acknowledgements.

Work at Caltech was supported by the Center for Molecular Magnetic Quantum
Materials, an Energy Frontier Research Center funded by the U.S. Department of
Energy, Office of Science, Basic Energy Sciences under award no. DE-SC0019330.
The computations were conducted at the Resnick High Performance Computing
Center, a facility supported by the Resnick Sustainability Institute at the
California Institute of Technology. Work at Lund University was supported by
grants from the Swedish research council (projects 2018-05003 and 2022-04978).
The computations were performed on computer resources provided by the Swedish
National Infrastructure for Computing (SNIC) at Lunarc at Lund University and
HPC2N at Umeå University, partially funded by the Swedish Research Council
(grant 2018-05973).

## Conflict of Interest

The authors have no conflicts to disclose.

## Data Availability

The data presented in this work can be reproduced using the open-source code
PySCF 2.0.1[46, 47] and block2 0.5.1.[37] The reference input and output files
can be found in the GitHub repository at https://github.com/hczhai/fe-dimer-
data.

## References

  * Beinert et al. [1997] Beinert, H.; Holm, R. H.; Munck, E. Iron-sulfur clusters: nature’s modular, multipurpose structures. _Science_ 1997, _277_ , 653–659.
  * Tanifuji and Ohki [2020] Tanifuji, K.; Ohki, Y. Metal–sulfur compounds in N2 reduction and nitrogenase-related chemistry. _Chemical reviews_ 2020, _120_ , 5194–5251.
  * Seefeldt et al. [2020] Seefeldt, L. C.; Yang, Z.-Y.; Lukoyanov, D. A.; Harris, D. F.; Dean, D. R.; Raugei, S.; Hoffman, B. M. Reduction of substrates by nitrogenases. _Chemical reviews_ 2020, _120_ , 5082–5106.
  * Rutledge and Tezcan [2020] Rutledge, H. L.; Tezcan, F. A. Electron transfer in nitrogenase. _Chemical reviews_ 2020, _120_ , 5158–5193.
  * Kirn and Rees [1992] Kirn, J.; Rees, D. Crystallographic structure and functional implications of the nitrogenase molybdenum–iron protein from Azotobacter vinelandii. _Nature_ 1992, _360_ , 553–560.
  * Peters et al. [1997] Peters, J. W.; Stowell, M. H.; Soltis, S. M.; Finnegan, M. G.; Johnson, M. K.; Rees, D. C. Redox-dependent structural changes in the nitrogenase P-cluster. _Biochemistry_ 1997, _36_ , 1181–1187.
  * Einsle et al. [2002] Einsle, O.; Tezcan, F. A.; Andrade, S. L.; Schmid, B.; Yoshida, M.; Howard, J. B.; Rees, D. C. Nitrogenase MoFe-protein at 1.16 Å resolution: a central ligand in the FeMo-cofactor. _Science_ 2002, _297_ , 1696–1700.
  * Spatzal et al. [2011] Spatzal, T.; Aksoyoglu, M.; Zhang, L.; Andrade, S. L.; Schleicher, E.; Weber, S.; Rees, D. C.; Einsle, O. Evidence for interstitial carbon in nitrogenase FeMo cofactor. _Science_ 2011, _334_ , 940–940.
  * Lancaster et al. [2011] Lancaster, K. M.; Roemelt, M.; Ettenhuber, P.; Hu, Y.; Ribbe, M. W.; Neese, F.; Bergmann, U.; DeBeer, S. X-ray emission spectroscopy evidences a central carbon in the nitrogenase iron-molybdenum cofactor. _Science_ 2011, _334_ , 974–977.
  * Jafari et al. [2022] Jafari, S.; Tavares Santos, Y. A.; Bergmann, J.; Irani, M.; Ryde, U. Benchmark Study of Redox Potential Calculations for Iron–Sulfur Clusters in Proteins. _Inorganic chemistry_ 2022, _61_ , 5991–6007.
  * Benediktsson and Bjornsson [2022] Benediktsson, B.; Bjornsson, R. Analysis of the Geometric and Electronic Structure of Spin-Coupled Iron–Sulfur Dimers with Broken-Symmetry DFT: Implications for FeMoco. _Journal of chemical theory and computation_ 2022, _18_ , 1437–1457.
  * Hoeke et al. [2019] Hoeke, V.; Tociu, L.; Case, D. A.; Seefeldt, L. C.; Raugei, S.; Hoffman, B. M. High-resolution ENDOR spectroscopy combined with quantum chemical calculations reveals the structure of nitrogenase Janus intermediate E4 (4H). _Journal of the American Chemical Society_ 2019, _141_ , 11984–11996.
  * Raugei et al. [2018] Raugei, S.; Seefeldt, L. C.; Hoffman, B. M. Critical computational analysis illuminates the reductive-elimination mechanism that activates nitrogenase for N 2 reduction. _Proceedings of the National Academy of Sciences_ 2018, _115_ , E10521–E10530.
  * Rohde et al. [2018] Rohde, M.; Sippel, D.; Trncik, C.; Andrade, S. L. A.; Einsle, O. The Critical E4 State of Nitrogenase Catalysis. _Biochemistry_ 2018, _57_ , 5497–5504, doi: 10.1021/acs.biochem.8b00509.
  * Siegbahn [2016] Siegbahn, P. E. M. Model Calculations Suggest that the Central Carbon in the FeMo-Cofactor of Nitrogenase Becomes Protonated in the Process of Nitrogen Fixation. _Journal of the American Chemical Society_ 2016, _138_ , 10485–10495.
  * Dance [2020] Dance, I. Computational Investigations of the Chemical Mechanism of the Enzyme Nitrogenase. _ChemBioChem_ 2020, _21_ , 1671–1709.
  * Cao et al. [2018] Cao, L.; Caldararu, O.; Ryde, U. Protonation and reduction of the FeMo cluster in nitrogenase studied by quantum mechanics/molecular mechanics (QM/MM) calculations. _Journal of chemical theory and computation_ 2018, _14_ , 6653–6678.
  * Cao and Ryde [2020] Cao, L.; Ryde, U. What Is the Structure of the E4 Intermediate in Nitrogenase? _Journal of chemical theory and computation_ 2020, _16_ , 1936–1952.
  * Thorhallsson et al. [2019] Thorhallsson, A. T.; Benediktsson, B.; Bjornsson, R. A model for dinitrogen binding in the E 4 state of nitrogenase. _Chemical science_ 2019, _10_ , 11110–11124.
  * Lovell et al. [2001] Lovell, T.; Li, J.; Liu, T.; Case, D. A.; Noodleman, L. FeMo Cofactor of Nitrogenase: A Density Functional Study of States $\mathrm{M^{N}}$, $\mathrm{M^{OX}}$, $\mathrm{M^{R}}$, and $\mathrm{M^{I}}$. _Journal of American Chemical Society_ 2001, _123_ , 12392–12410.
  * Cao and Ryde [2018] Cao, L.; Ryde, U. Influence of the protein and DFT method on the broken-symmetry and spin states in nitrogenase. _International Journal of Quantum Chemistry_ 2018, _118_ , e25627.
  * Hohenberg and Kohn [1964] Hohenberg, P.; Kohn, W. Inhomogeneous electron gas. _Physical Review_ 1964, _136_ , B864.
  * Cao and Ryde [2019] Cao, L.; Ryde, U. Extremely large differences in DFT energies for nitrogenase models. _Physical Chemistry Chemical Physics_ 2019, _21_ , 2480–2488.
  * Wei and Siegbahn [2021] Wei, W.-J.; Siegbahn, P. E. The active E4 structure of nitrogenase studied with different DFT functionals. _Journal of computational chemistry_ 2021, _42_ , 81–85.
  * White [1992] White, S. R. Density matrix formulation for quantum renormalization groups. _Physical review letters_ 1992, _69_ , 2863.
  * White [1993] White, S. R. Density-matrix algorithms for quantum renormalization groups. _Physical Review B_ 1993, _48_ , 10345.
  * Chan and Head-Gordon [2002] Chan, G. K.-L.; Head-Gordon, M. Highly correlated calculations with a polynomial cost algorithm: A study of the density matrix renormalization group. _The Journal of chemical physics_ 2002, _116_ , 4462–4476.
  * Chan and Sharma [2011] Chan, G. K.-L.; Sharma, S. The density matrix renormalization group in quantum chemistry. _Annual review of physical chemistry_ 2011, _62_ , 465–481.
  * Mitrushenkov et al. [2001] Mitrushenkov, A. O.; Fano, G.; Ortolani, F.; Linguerri, R.; Palmieri, P. Quantum chemistry using the density matrix renormalization group. _The Journal of Chemical Physics_ 2001, _115_ , 6815–6821.
  * Chan et al. [2016] Chan, G. K.-L.; Keselman, A.; Nakatani, N.; Li, Z.; White, S. R. Matrix product operators, matrix product states, and ab initio density matrix renormalization group algorithms. _The Journal of chemical physics_ 2016, _145_ , 014102.
  * Sharma and Chan [2012] Sharma, S.; Chan, G. K.-L. Spin-adapted density matrix renormalization group algorithms for quantum chemistry. _The Journal of chemical physics_ 2012, _136_ , 124121.
  * Olivares-Amaya et al. [2015] Olivares-Amaya, R.; Hu, W.; Nakatani, N.; Sharma, S.; Yang, J.; Chan, G. K.-L. The ab-initio density matrix renormalization group in practice. _The Journal of chemical physics_ 2015, _142_ , 034102.
  * Wouters and Van Neck [2014] Wouters, S.; Van Neck, D. The density matrix renormalization group for ab initio quantum chemistry. _The European Physical Journal D_ 2014, _68_ , 272.
  * Keller et al. [2015] Keller, S.; Dolfi, M.; Troyer, M.; Reiher, M. An efficient matrix product operator representation of the quantum chemical Hamiltonian. _The Journal of chemical physics_ 2015, _143_ , 244118.
  * Baiardi and Reiher [2020] Baiardi, A.; Reiher, M. The density matrix renormalization group in chemistry and molecular physics: Recent developments and new challenges. _The Journal of Chemical Physics_ 2020, _152_ , 040903.
  * Brabec et al. [2021] Brabec, J.; Brandejs, J.; Kowalski, K.; Xantheas, S.; Legeza, Ö.; Veis, L. Massively parallel quantum chemical density matrix renormalization group method. _Journal of Computational Chemistry_ 2021, _42_ , 534–544.
  * Zhai and Chan [2021] Zhai, H.; Chan, G. K.-L. Low communication high performance ab initio density matrix renormalization group algorithms. _The Journal of Chemical Physics_ 2021, _154_ , 224116.
  * Sharma et al. [2014] Sharma, S.; Sivalingam, K.; Neese, F.; Chan, G. K. Low-energy spectrum of iron–sulfur clusters directly from many-particle quantum mechanics. _Nature chemistry_ 2014, _6_ , 927–933.
  * Li et al. [2019] Li, Z.; Guo, S.; Sun, Q.; Chan, G. K. Electronic landscape of the P-cluster of nitrogenase as revealed through many-electron quantum wavefunction simulations. _Nature chemistry_ 2019, _11_ , 1026–1033.
  * Li et al. [2019] Li, Z.; Li, J.; Dattani, N. S.; Umrigar, C.; Chan, G. K.-L. The electronic complexity of the ground-state of the FeMo cofactor of nitrogenase as relevant to quantum simulations. _The Journal of chemical physics_ 2019, _150_ , 024302.
  * Tao et al. [2003] Tao, J.; Perdew, J. P.; Staroverov, V. N.; Scuseria, G. E. Climbing the density functional ladder: Nonempirical meta–generalized gradient approximation designed for molecules and solids. _Physical Review Letters_ 2003, _91_ , 146401.
  * Schäfer et al. [1992] Schäfer, A.; Horn, H.; Ahlrichs, R. Fully optimized contracted Gaussian basis sets for atoms Li to Kr. _The Journal of Chemical Physics_ 1992, _97_ , 2571–2577.
  * Grimme et al. [2010] Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu. _The Journal of chemical physics_ 2010, _132_ , 154104.
  * Bartlett and Musiał [2007] Bartlett, R. J.; Musiał, M. Coupled-cluster theory in quantum chemistry. _Reviews of Modern Physics_ 2007, _79_ , 291.
  * Shavitt and Bartlett [2009] Shavitt, I.; Bartlett, R. J. _Many-body methods in chemistry and physics: MBPT and coupled-cluster theory_ ; Cambridge university press, 2009.
  * Sun et al. [2018] Sun, Q.; Berkelbach, T. C.; Blunt, N. S.; Booth, G. H.; Guo, S.; Li, Z.; Liu, J.; McClain, J. D.; Sayfutyarova, E. R.; Sharma, S., et al. PySCF: the Python-based simulations of chemistry framework. _Wiley Interdisciplinary Reviews: Computational Molecular Science_ 2018, _8_ , e1340.
  * Sun et al. [2020] Sun, Q.; Zhang, X.; Banerjee, S.; Bao, P.; Barbry, M.; Blunt, N. S.; Bogdanov, N. A.; Booth, G. H.; Chen, J.; Cui, Z.-H., et al. Recent developments in the PySCF program package. _The Journal of chemical physics_ 2020, _153_ , 024109.
  * Saue [2011] Saue, T. Relativistic Hamiltonians for chemistry: A primer. _ChemPhysChem_ 2011, _12_ , 3077–3094.
  * Peng and Reiher [2012] Peng, D.; Reiher, M. Exact decoupling of the relativistic Fock operator. _Theoretical Chemistry Accounts_ 2012, _131_ , 1–20.
  * Kutzelnigg and Liu [2005] Kutzelnigg, W.; Liu, W. Quasirelativistic theory equivalent to fully relativistic theory. _The Journal of chemical physics_ 2005, _123_ , 241102.
  * Becke [1988] Becke, A. D. Density-functional exchange-energy approximation with correct asymptotic behavior. _Physical review A_ 1988, _38_ , 3098\.
  * Lee et al. [1988] Lee, C.; Yang, W.; Parr, R. G. Development of the Colle-Salvetti correlation-energy formula into a functional of the electron density. _Physical review B_ 1988, _37_ , 785.
  * Becke [1993] Becke, A. D. A new mixing of Hartree–Fock and local density-functional theories. _The Journal of chemical physics_ 1993, _98_ , 1372–1377.
  * Schurkus et al. [2020] Schurkus, H.; Chen, D.-T.; Cheng, H.-P.; Chan, G.; Stanton, J. Theoretical prediction of magnetic exchange coupling constants from broken-symmetry coupled cluster calculations. _The Journal of Chemical Physics_ 2020, _152_ , 234115.
  * Yamaguchi et al. [1986] Yamaguchi, K.; Takahara, Y.; Fueno, T. Ab-initio molecular orbital studies of structure and reactivity of transition metal-oxo compounds. Applied Quantum Chemistry: Proceedings of the Nobel Laureate Symposium on Applied Quantum Chemistry in Honor of G. Herzberg, RS Mulliken, K. Fukui, W. Lipscomb, and R. Hoffman, Honolulu, HI, 16–21 December 1984. 1986; pp 155–184.
  * Pipek and Mezey [1989] Pipek, J.; Mezey, P. G. A fast intrinsic localization procedure applicable for abinitio and semiempirical linear combination of atomic orbital wave functions. _The Journal of Chemical Physics_ 1989, _90_ , 4916–4926.
  * Lee et al. [2021] Lee, S.; Zhai, H.; Sharma, S.; Umrigar, C. J.; Chan, G. K.-L. Externally Corrected CCSD with Renormalized Perturbative Triples (R-ecCCSD (T)) and the Density Matrix Renormalization Group and Selected Configuration Interaction External Sources. _Journal of Chemical Theory and Computation_ 2021, _17_ , 3414–3425.
  * Møller and Plesset [1934] Møller, C.; Plesset, M. S. Note on an approximation treatment for many-electron systems. _Physical review_ 1934, _46_ , 618\.
  * Pople et al. [1976] Pople, J. A.; Binkley, J. S.; Seeger, R. Theoretical models incorporating electron correlation. _International Journal of Quantum Chemistry_ 1976, _10_ , 1–19.
  * Neese and Valeev [2011] Neese, F.; Valeev, E. F. Revisiting the atomic natural orbital approach for basis sets: robust systematic basis sets for explicitly correlated and conventional correlated ab initio methods? _Journal of chemical theory and computation_ 2011, _7_ , 33–43.
  * Perdew et al. [1996] Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized gradient approximation made simple. _Physical review letters_ 1996, _77_ , 3865.
  * Grimme [2006] Grimme, S. Semiempirical GGA-type density functional constructed with a long-range dispersion correction. _Journal of computational chemistry_ 2006, _27_ , 1787–1799.
  * Furness et al. [2020] Furness, J. W.; Kaplan, A. D.; Ning, J.; Perdew, J. P.; Sun, J. Accurate and numerically efficient r2SCAN meta-generalized gradient approximation. _The journal of physical chemistry letters_ 2020, _11_ , 8208–8215.
  * Staroverov et al. [2003] Staroverov, V. N.; Scuseria, G. E.; Tao, J.; Perdew, J. P. Comparative assessment of a new nonempirical density functional: Molecules and hydrogen-bonded complexes. _The Journal of chemical physics_ 2003, _119_ , 12129–12137.
  * Salomon et al. [2002] Salomon, O.; Reiher, M.; Hess, B. A. Assertion and validation of the performance of the B3LYP* functional for the first transition metal row and the G2 test set. _The Journal of chemical physics_ 2002, _117_ , 4729–4737.
  * Perdew et al. [1996] Perdew, J. P.; Ernzerhof, M.; Burke, K. Rationale for mixing exact exchange with density functional approximations. _The Journal of chemical physics_ 1996, _105_ , 9982–9985.
  * Zhao and Truhlar [2008] Zhao, Y.; Truhlar, D. G. The M06 suite of density functionals for main group thermochemistry, thermochemical kinetics, noncovalent interactions, excited states, and transition elements: two new functionals and systematic testing of four M06-class functionals and 12 other functionals. _Theoretical chemistry accounts_ 2008, _120_ , 215–241.
  * Cancès et al. [2013] Cancès, E.; Maday, Y.; Stamm, B. Domain decomposition for implicit solvation models. _The Journal of Chemical Physics_ 2013, _139_ , 054111.

## TOC Graphic