We first examine the effect of including the relativistic correction in the framework of Douglas-Kroll-Hess (DKH) to the DFT-D3(BJ) interaction energy and showed the results in Table . In general, the DKH method lowers the interaction energy by no more than 0.01 eV. As the improvement is very small, further calculations will not employ relativistic treatment. Next, we investigated the effect of using different sizes of DFT integration less than a half of the number of basis functions than the former, we proceed our interaction energy study with the latter (ma-def2-TZVP) basis set. In Table , we show that the self-consistent and non self-consistent computation of the non-local dispersion correction yield essentially the same results (within 0.003 eV). Therefore, the results of the computationally inexpensive non self-consistent DFT-NL computation is the ones shown in the manuscript. Table shows DLPNO-CCSD(T) interaction energy using def2 and ma-def2 series of basis set. We first notice that the basis-set superposition error is larger in ma-def2-QZVP than in def2-QZVP calculations, and hence, the results obtained from the def2-QZVP calculations serve as the reference. Secondly, the cheaper NormalPNO settings of the DLPNO-CCSD(T) method result in slightly lower def2-QZVP interaction energy than that of the TightPNO settings by 0.009 eV. Finally, the def2-TZVP calculations slightly overbind by 0.017 eV, as compared to the calculations using the def2-QZVP basis set. Since a combination of loosening the PNO threshold from TightPNO to NormalPNO and reducing the basis set size from def2-QZVP to def2-TZVP presents a combination of parameters that only give a binding energy that is higher by 0.008 eV from the reference, we will employ this recipe for computing the DLPNO-CCSD(T) interaction energy of the larger [Fe(bpy