The calculations of the adiabatic energy difference and the binding energies are performed using Quantum Espresso. We employ the GBRV ultrasoft pseudopotentials and wavefunction and charge density cutoffs of 70 Ry and 700 Ry, respectively. The reciprocal space sampling is performed by using a 2×2×1 Monkhorst-Pack grid. A full geometrical optimization is performed for both the bare and the loaded MOF using PBE+U with U=4 on both the octahedrally coordinated Fe and the open metal center M (vide infra), starting from the published structure. In order to compute the spin state energetics we employ instead a Hubbard U computed self-consistently 45 using linear-response, called here U sc. The U sc is computed separately for each metal and spin state by employing the optimized geometry. For Fe and Ni we employ a molecular fragment M(CNH)4 and two N2 adsorbed molecules. For Mn, the U sc is computed using the periodic Fe(pz)[Mn(CN)4] and one CO molecule. We employ ortho-atomic projectors. These U sc are then employed for all the subsequent calculations including those for the bare MOF and the other gas molecules, for both Fe(bpac)M(CN)4 and Fe(pz)M(CN)4. We note several convergence issues in the self-consistent cycle when using U sc=8.9 eV for the octahedrally-coordinated Fe. Thus, in order to avoid convergence problems we have set U =4 eV for this atom for every MOF calculation. Similar to the molecular case, all the PBE+U calculations are performed using the Grimme-D3 correction together with the BJ damping scheme. As reported in previous works, the Hubbard correction is needed in order to correctly describe the transition-metal 3d electrons in MOFs.