These findings highlight that the periodic models with both types of axial ligands predict an HS configuration for both pyridinic and pyrrolic structures when including proper correlation effects in the calculations. Each relaxed structure shows substantial displacement of iron from the carbon layer, similar to the OOP phase of the square-planar models, which stabilizes the calculated spin configurations. In contrast, the more standard approach of using a pure GGA formalism consistently prefers an IS state associated with smaller iron displacements. This emphasizes that the structural composition of the active site model as well as the computational approach need to be developed carefully for meaningful comparisons with experimental data.