Testing Multireference Methods Against a Measured Iron Spin Gap: the Peters Fe₂(μ-H)₂ Complex

Abstract

Spin States of a Diiron Dihydride — a Computational Case Study of the Peters 3-N₂ Complex In plain language We asked computer simulation to reproduce one measured property of a real iron compound, and got the direction right but the size only roughly. The compound is a two-iron molecule made at Caltech in 2014. Its electrons can sit fully paired, or with two unpaired, and the lab measured the paired arrangement as lower in energy by 4.7 kcal/mol — a fraction of a chemical bond. Why it matters. Nature fixes nitrogen from the air using an iron-rich cluster in an enzyme; industry needs the energy-hungry Haber–Bosch process. Understanding the enzyme starts with simulations that can be trusted on iron, where simple methods are often wrong by tens of kcal/mol. So every method must first pass tests like this one, where the answer is known. What we found. Our careful method picks the right winner: the paired arrangement is lower. On the size of the gap there is a subtlety. Each arrangement prefers a slightly different shape of the molecule, and the fair comparison is each in its own best shape — like comparing two people's heights with both standing. Measured in one shape, our gap is 8.2 kcal/mol; in the other, 11.4. Allowing each its own shape, our best estimate is about +1, against the measured 4.7. The puzzle. When the molecule changes shape, our careful method and a simpler one disagree about which arrangement benefits, by about 20 kcal/mol. We tested and ruled out the obvious causes one by one: the simplified molecule, the handling of the shapes, the level of detail, the innermost electrons, and the weak sticky forces between atoms. We also ruled out being lost in the fog. These calculations find the lowest energy a bit like walking downhill in fog, where you can stop in a side valley and believe you have reached the bottom while still high up. Restarting the walk from a completely different place led to exactly the same spot. What it probably is, and what it needs. The disagreement comes from the methods themselves. The likeliest culprits are two known blind spots. To stay affordable, our careful method treats only a small group of electrons in full detail, and iron is known to need more. It also handles iron's inner-shell electrons only approximately, which is known to favour the unpaired arrangement — exactly the way our answer misses. Settling it needs heavier computing than a home desktop: either a much larger fully treated group, or an independent reference method run on data-centre hardware. We tried a smaller version of that reference test at home; it was too slow and only borderline trustworthy for this molecule. What else came out of it. Along the way we built a pipeline that splits work between a graphics card and the main processor, checks every handover, and refuses to continue when a result fails a sanity check. Those checks stopped several real errors — a wrong starting solution, distorted orbitals, mismatched comparisons — before they could reach a result. Abstract State-specific CASSCF/SC-NEVPT2 in a balanced CAS(8e, 8o) reproduces the singlet ground state of the Peters diiron dihydride 3-N₂ but not yet the size of its singlet–triplet gap. The experimental gap is 4.7 kcal/mol (19.8 kJ/mol), from a two-state fit to the variable-temperature ¹H NMR shift of the bridging hydrides (Rittle, McCrory and Peters, 2014). On the published BP86 geometries, NEVPT2 vertical gaps are +8.20 kcal/mol at the singlet geometry and +11.37 at the triplet geometry. A composite adiabatic gap — NEVPT2 vertical plus BP86 geometric relaxation, by two routes — is +0.6 to +1.7 kcal/mol (mean +1.1). The methods disagree on how the gap changes between the two geometries: BP86 has it falling by 17.3 kcal/mol, CASSCF and NEVPT2 rising by about 3. We excluded, by direct test, the ligand trimming and geometry handling (BP86 on the model reproduces the full ligand within 0.5 kcal/mol), the basis and core correlation (MP2 changes by under 1 kcal/mol across four treatments), NEVPT2-specific artefacts (MP2 reproduces NEVPT2's geometry dependence for the singlet), dispersion (D3 has the opposite sign), and local minima (an independent restart reached the identical CASSCF energy). Two documented candidates remain. The active space lacks a correlating 3d′ shell, and second-order perturbation theory describes metal 3s3p semicore correlation erratically, a known bias toward high spin that matches the direction of our 3.6 kcal/mol shortfall. Testing either needs a much larger active space or coupled-cluster on the full model, beyond a home workstation. All calculations used PySCF 2.14 with GPU4PySCF 1.8.1 on an RTX 5090. Every CPU–GPU handover was verified numerically (agreement ~10⁻¹⁰ Ha), and run-time invariant checks refused several silently wrong intermediate results. Background and the experimental target The target is the singlet–triplet gap of 3-N₂, measured at +4.7 kcal/mol with the singlet lower. 3-N₂ is the Fe(II)Fe(II) member of a family of Fe₂(μ-H)₂ complexes whose N₂ affinity rises a million-fold on one-electron reduction (Rittle, McCrory and Peters, J. Am. Chem. Soc. 2014, 136, 13853). Two bridging hydrides join the irons. One dinucleating ligand spans both: each iron binds one silyl silicon and two phosphines, and the two silicons are linked through an oxygen. One N₂ is bound. The measurement. The bridging-hydride ¹H NMR shift varies strongly with temperature. The supporting information fits it to a singlet ground state with a thermally accessible triplet, giving an energy difference of 19,800 ± 94 J/mol (4.73 kcal/mol). The quoted error is the fit's alone; with the two-state model's own assumptions, a realistic uncertainty is about ±1 kcal/mol. A fit of this kind measures an energy difference between thermally relaxed states, so the comparison is to an adiabatic gap. Why it is hard. Spin-state gaps of first-row transition-metal complexes are a known weak point of electronic-structure methods. Density functionals can disagree by tens of kcal/mol, and Hartree–Fock here puts the triplet below the singlet. Multireference methods treat the near-degeneracy properly, but only within a chosen active space, and the dynamic correlation outside it must be added perturbatively. Geometries. The supporting information gives BP86/6-31G(d) geometries, with the full ligand, for 3-N₂ as singlet and triplet. Where experiment exists for the family, the published DFT reproduces the Fe–Fe distance within 0.03 Å: 2.494 against 2.480 Å (X-ray, 3-N₂), 2.729 against 2.755 Å (EXAFS, 3-(N₂)₂), 2.839 against 2.872 Å (X-ray, 4-(N₂)₂). The two 3-N₂ geometries differ by 0.07 Å in Fe–Fe, and, after rigid alignment, by 0.30 Å RMSD over the whole molecule: ligand heavy atoms move by up to 0.83 Å. Models, geometries and methods The production model is the published full-ligand structure with its eight P-isopropyl groups cut to methyl: 85 atoms, def2-SVP on the 13 atoms that matter for the spin question, STO-3G elsewhere. Each stage below was validated before use; versions and hashes are in the data section. Geometries. The two published BP86/6-31G(d) structures of 3-N₂ (singlet and triplet), extracted from the supporting-information PDF. The triplet geometry is rigidly rotated onto the singlet's frame (Kabsch; RMSD 0.435 → 0.299 Å), so orbitals can be carried between geometries. Energies are unaffected. Model. Each isopropyl on phosphorus is cut to methyl, the removed carbons replaced by hydrogens along the original bonds at 1.089 Å. A 29-atom "core" model — iron, hydrides, N₂ and the P/Si/O donors, hydrogen-capped — was used only for tests. Basis. def2-SVP on Fe, the hydrides, N, P, Si and O; STO-3G on the 72 outer C and H: 422 basis functions instead of 870, about 9× cheaper per slow-path integral build. Density fitting, with one shared integral set per geometry. SCF. RHF (singlet) and ROHF (triplet), run on the GPU with GPU4PySCF, then re-evaluated on the CPU. Fallback: CPU DIIS and Newton. An internal stability check is applied once. At the triplet geometry each spin starts from its own converged orbitals at the singlet geometry, projected. Active space. AVAS projection onto Fe 3d and hydride 1s functions, at a fixed size of 12 orbitals, gives CAS(16e, 12o) for both spins. A 60-cycle 12-orbital singlet CASSCF (oscillating, not converged) left four orbitals doubly occupied to three decimals in both spins (subspace overlap 0.9997). These move to the core, giving CAS(8e, 8o). At the triplet geometry the space is carried over by AVAS with the converged orbitals as target (match 0.986 chosen, 0.014 left out). Correlated energies. State-specific CASSCF(8e, 8o), held to a pure spin by a penalty (⟨S²⟩ 0.0000 and 2.0000), then strongly contracted NEVPT2, density-fitted, all electrons correlated. On O₂ the same code gives a singlet–triplet gap of 23.35 kcal/mol (experiment 22.6). Composite. BP86/6-31G(d) single points (RKS/UKS) on the full published structures supply the geometric relaxation energies. Diagnostics. A CAS(12e, 12o) CASCI probe, DF-MP2 in four treatments, D3(BJ) dispersion, and a same-state restart; described with their results below. Hardware. One workstation: Intel Core Ultra 9 285K (24 threads), NVIDIA RTX 5090 (32 GB), under WSL2. The GPU ran the SCF and the integral builds inside CASSCF, about 4.5 s per build against 9.5 s on the CPU. CI solves, NEVPT2 and stability analysis ran on the CPU. Each session was held to a 1.8-hour budget, with checkpoints. Results: the singlet–triplet gap NEVPT2 puts the singlet below the triplet at both geometries, as experiment does; the composite adiabatic gap, +1.1 kcal/mol, falls 3.6 short of the measured +4.7. Figure 1 — The composite adiabatic gap is +1.1 kcal/mol; experiment is +4.7. Singlet–triplet gap by method: vertical gaps at both geometries, adiabatic gaps, and the experimental band (file figures/fig1_gap_by_method.png in this record). Hartree–Fock gets the order wrong a