of 11
Multireference
Protonation
Energetics
of a Dimeric
Model
of
Nitrogenase
Iron
Sulfur
Clusters
Published
as part of The Journal
of Physical
Chemistry
A virtual
special
issue “Roland
Lindh
Festschrift”.
Huanchen
Zhai,
*
Seunghoon
Lee, Zhi-Hao
Cui, Lili Cao,
Ulf Ryde,
*
and Garnet
Kin-Lic
Chan
*
Cite This:
J. Phys.
Chem.
A
2023,
127, 9974−9984
Read Online
ACCESS
Metrics
& More
Article
Recommendations
*
Supporting
Information
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
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.
1. INTRODUCTION
Nitrogenase
is the only enzyme
that can catalyze
the
conversion
of atmospheric
dinitrogen
(N
2
) to ammonia
(NH
3
), the key reaction
in nitrogen
fixation.
1
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
], which
transfers
electrons
to the FeMoco
active site.
5,6
During
the past two
decades,
the atomic
structures
of the nitrogenase
clusters
have
been determined
by X-ray crystallography.
7
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
E
4
intermediate
of nitrogenase,
12
formed
by adding
four electrons
and protons
to the E
0
resting
state of FeMoco
and responsible
for the binding
of N
2
, has
been studied
computationally
by several
groups.
13
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 E
4
state).
11,23,24
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
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 nonclassical
spin correla-
tions.
38
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
Received:
September
12, 2023
Revised:
October
31, 2023
Accepted:
November
1, 2023
Published:
November
15,
2023
Article
pubs.acs.org/JPCA
© 2023
The Authors.
Published
by
American
Chemical
Society
9974
https://doi.org/10.1021/acs.jpca.3c06142
J. Phys.
Chem.
A
2023,
127, 9974
9984
This article is licensed under CC-BY 4.0
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.
2. 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.
17
We denote
the four
protonated
configurations
HC, HS, HFe, and HFe
2
,
respectively.
Our goal of the study is to identify
the lowest-
energy
structure
and to predict
the relative
energetics
of the
four structures.
2.1. 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 performed
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
Figure
1 and give the coordinates
in the Supporting
Information
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 1. 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 length
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
geo-
Figure
1.
Geometries
of the protonated
dimer complex
[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.
Table
1. Bond
Lengths
(in Å) between
the
Iron
Ions
(Fe1
and
Fe2)
and
Their
Direct
Ligands
in the
Optimized
Geometries
of the Four
Dimer
Models
a
bond
HC
HS
HFe
HFe
2
S
t1
Fe1
2.45
2.43
2.62
2.43
S
t2
Fe1
2.45
2.45
2.50
2.48
S
t3
Fe2
2.45
2.44
2.42
2.39
S
t4
Fe2
2.42
2.39
2.38
2.43
S
br
Fe1
2.26
2.47
2.42
2.31
S
br
Fe2
2.28
2.53
2.21
2.34
average
(all S
Fe)
2.38
2.45
2.42
2.40
C
br
Fe1
2.20
1.98
2.02
1.97
C
br
Fe2
2.19
1.99
1.95
1.99
H
Fe1
1.76
1.72
H
Fe2
1.76
a
S
t1
S
t4
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
2
ion,
and H is the added proton,
when binding
to Fe.
The Journal
of Physical
Chemistry
A
pubs.acs.org/JPCA
Article
https://doi.org/10.1021/acs.jpca.3c06142
J. Phys.
Chem.
A
2023,
127, 9974
9984
9975
metric 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 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.
2.2. 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
E
E
E
E
(
)
composite
CCSD(T)
TZ
DM
RG
act
CC
SD(T)
act
CCSD(T)
CB
S
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 multi-
reference
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,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 performed
additional
calcu-
lations
to understand
the impact
of different
approximations
and to estimate
the reliability
of the corrections.
We describe
the different
terms and these aspects
below.
2.3. Coupled
Cluster
Calculations.
For the coupled
cluster (CC) calculations,
we started
from the BS unrestricted
reference
determinants.
To understand
the influence
of orbitals
and importance
of triples, we first performed
calculations
using
the small cc-pVDZ-DK
basis set and the exact two-component
(X2C)
scalar relativistic
Hamiltonian
48
50
(larger
basis set CC
calculations
are discussed
in the basis set correction
section
below)
and using 40 frozen core orbitals
to reduce
computa-
tional cost. To examine
the impact
of the orbital
choice,
we
used both Kohn
Sham
DFT (with the TPSS or B3LYP
functionals
51
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
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,45
For the DFT reference
determinants,
we computed
CCSD(T)
results
based on the semi-canon-
icalized
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
=
=
J
E
E
S
S
E
E
S
S
HS
BS
BS
2
HS
2
PS
BS
BS
2
PS
2
where
J
is the exchange
coupling.
For simplicity,
we used
S
PS
2
= 0 and
S
BS
2
and
S
HS
2
computed
at the CCSD
level
in the preceding
formula
for computing
J
from both CCSD and
CCSD(T)
energies.
The difference
between
the BS and PS
states 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.
2.4. 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 performed
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
the Supporting
Information
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 full space into
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 on the cc-pVDZ-DK
basis
contains
180 electrons
in 321 spatial orbitals.
We constructed
the active spaces from the UHF/CCSD
natural
orbitals
(see
the Supporting
Information
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
of 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 Journal
of Physical
Chemistry
A
pubs.acs.org/JPCA
Article
https://doi.org/10.1021/acs.jpca.3c06142
J. Phys.
Chem.
A
2023,
127, 9974
9984
9976
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 Table 2.
2.5. Basis Set Correction
and Relativistic
Contribu-
tion.
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
=
E
X
E
Y
E
X
Y
X
Y
corr
corr
(
)
corr
(
)
where
X
= 2 (DZ) and
Y
= 3 (TZ) and taking
β
= 2.4. For the
corresponding
mean-field
energy
at the CBS limit, we simply
used
E
UHF
=
E
UHF
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
|
|
E
E
E
E
1
2
(
)
CCSD(T)
DZ / TZ
CCSD(T)
exact
M
P2
TZ / QZ
M
P2
DZ / TZ
where
Δ
E
CCSD(T)
DZ/TZ
→∞
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.
2.6. 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.
43
2.7. 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-decomposi-
tion COSMO
solvation
model
68
with a dielectric
constant
ε
=
4.0. We compared
the relative
energies
with and without
the
Table
2. Relative
Single-Point
Energies
(in kJ/mol)
for the
BS State
of the
Four
Protonated
Fe Dimer
Structures
Computed
by Using
Different
Theories
with
the
cc-PVDZ-
DK
Basis
and
the
Scalar
Relativistic
X2C
Hamiltonian
a
energy
difference
E
E
HC
theory
HC
HS
HFe
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
a
The energy
of the HC structure
is used as the reference.
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 structure
is used as the reference.
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
a
energy
difference
E
E
HC
theory
HC
HS
HFe
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
S
2
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
S
2
3.89
3.78
4.57
4.15
exchange
coupling
J
(estimated,
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
a
The energy
of the HC structure
is used as the reference.
The Journal
of Physical
Chemistry
A
pubs.acs.org/JPCA
Article
https://doi.org/10.1021/acs.jpca.3c06142
J. Phys.
Chem.
A
2023,
127, 9974
9984
9977
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 3, for the four structures.
Solvation
also leads to a modest,
but nonuniform,
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
HFe
d
2
E
HC
. Clearly
solvation
is
important
for accurate
studies
of 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.
3. RESULTS
AND DISCUSSION
In Section
3.1, we first discuss
the CC energies
on 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
3.2 we discuss
the contribution
associated
with correcting
the BS spin states.
In Section
3.3 we discuss
the detailed
multireference
corrections
entering
the composite
energy
formula
from the
DMRG
and CC calculations.
In Section
3.4 we discuss
CC
calculations
in larger basis sets, the CBS extrapolation
entering
the composite
energy,
and the size of the relativistic
effects.
We
report
the final composite
energies,
the prediction
of the
lowest energy
protonation
site, the relative
ordering,
and the
comparison
with DFT calculations
in Section
3.5.
3.1. CC Energies:
Importance
of Higher-Order
Correlations.
We show the energies
of the four protonated
structures
relative
to the HC structure
from calculations
with
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), (55o, 48e), and (63o, 64e) active spaces.
The energy
of the HC structure
is used as the
reference.
The trends in the (T) and DMRG
corrections
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
a
energy
correction
E
E
CCSD
theory
HC
HS
HFe
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
a
The CCSD
energy
is used as a reference
for the 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)”.
The Journal
of Physical
Chemistry
A
pubs.acs.org/JPCA
Article
https://doi.org/10.1021/acs.jpca.3c06142
J. Phys.
Chem.
A
2023,
127, 9974
9984
9978
HF, CCSD,
and CCSD(T)
methods
with the cc-pVDZ-DK
basis in Table 2 and Figure 2. All methods
found that the HC
structure
is the most stable. However,
CCSD
and CCSD(T)
predict
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
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 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 2). 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.
3.2. Spin-State
Corrections.
As the preceding
calcula-
tions 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 of 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.
3.3. Multireference
Effects.
To obtain a more complete
picture
of the multireference
effects,
we next considered
ab
initio
DMRG
energies.
In Figures
3a and 3b 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
analysis
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 (Figure
3b), 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 Figure
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 Figure 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 HFe
2
is 60%, 66%, and 63% of the
contribution
in the full orbital space, respectively.
The largest
Figure
5.
CCSD(T)
and DMRG
correlation
energies,
relative
to the
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.
Table
5. UHF,
MP2
Correlation,
CCSD
Correlation,
and
(T)
Correction
Energies
Computed
Using
Different
Basis
Sets
a
E
E
average
(kJ/mol)
basis
E
average
(hartree)
HC
HS
HFe
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
a
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.
The Journal
of Physical
Chemistry
A
pubs.acs.org/JPCA
Article
https://doi.org/10.1021/acs.jpca.3c06142
J. Phys.
Chem.
A
2023,
127, 9974
9984
9979
estimated
error for the extrapolated
DMRG
energy (3 kJ/mol)
illustrates
that the DMRG
energies
are almost
exactly
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 Figure 4.
Based on the data mentioned
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
the 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 multi-
reference
character
of the structure
increases.
In Figure
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
of using 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.
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 kJ/mol
shift is added for clarity.
Table
6. Scalar
Relativistic
Corrections
(in kJ/mol)
to the
Relative
Energies
Computed
Using
Different
Theories
on
the def2-SV(P)
Basis
a
energy
difference
E
E
HC
theory
HC
HS
HFe
HFe
2
relativistic:
Δ
E
X2C
Δ
E
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
a
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.
The Journal
of Physical
Chemistry
A
pubs.acs.org/JPCA
Article
https://doi.org/10.1021/acs.jpca.3c06142
J. Phys.
Chem.
A
2023,
127, 9974
9984
9980
3.4. Basis Set Correction
and Relativistic
Contribu-
tion.
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 larger basis
sets. The results are listed in Table 5 and plotted
in Figure 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 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
nonmonotonically
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, Figure 6) is dominated
by the
basis set trend of the CCSD
relative
energies
beyond
TZ (top
left panel, Figure 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.
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 HFe
2
[11
12
kJ/mol
at the CCSD(T)
level].
Relativistic
effects are clearly necessary
to describe
differential
bonding
to Fe.
3.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 Figure 7. Overall,
we find that
E
HC
<
E
HS
<
E
HFe
d
2
<
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
Table
7. Relative
Single-Point
Energies
(in kJ/mol)
for the
Four
Protonated
Fe Dimer
Structures
Computed
by the
Composite
Method
a
energy
difference
E
E
HC
UKS deviation
correction/
functional
HC
HS
HFe
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
r
2
SCAN
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
a
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 maximum
deviation
in the
energy difference
E
E
HC
between
the composite
and UKS methods
for each functional.
Figure
7.
Comparison
between
the difference
in the 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.
The Journal
of Physical
Chemistry
A
pubs.acs.org/JPCA
Article
https://doi.org/10.1021/acs.jpca.3c06142
J. Phys.
Chem.
A
2023,
127, 9974
9984
9981