To identify whether a bent geometry is feasible in the periodic models, we manually placed the iron atom 0.85 Å above the carbon layer of the relaxed pyridinic and pyrrolic structures, and then fully relaxed them using GGA+U. A bent geometry is found only for the pyrrolic cases, while the pyridinic case prefers to remain flat. In this 'out-of-plane' (OOP) phase of the pyridinic model, the iron atom is found 1.15 Å above the layer, accompanied by a 7.3% increase in Fe-N distances to 1.99 Å compared to the flat 6×6 pyrrolic structures, see Figure . Within GGA+U, the most stable spin state in the 6×6 OOP graphene environment is a 2S + 1 = 5 high spin configuration: the d z 2 orbital is now doubly occupied, whereas the other orbitals remain singly occupied. According to our non-spinpolarized GGA results, unlike the flat 6×6 pyrrolic structure where a sizable energy gap of approximately 3eV is calculated between the highest-lying orbital (d xy ) and the second highestlying one (d yz ), this gap is significantly reduced to about 2 eV in the OOP phase, likely enabling a high-spin state. The vertical distance between the graphene layer and the iron atom is slightly reduced in the 4×4 case to 0.72 Å, whereas it increases to 1.46 Å in the 8×8 case. According to the GGA+U results, the OOP phase is energetically destabilised by 81 meV (≈ 1.87 kcal/mol) compared to the flat model. This difference is much larger (464 meV ≈ 10.7 kcal/mol) using GGA. Using the hybrid functional PBE0, a similar high spin configuration can be found with a different shape on the d-orbital of Fe atom (see Figure ). A detailed observation when comparing the electronic structure of pure-GGA PBE functional and PBE0 hybrid functionals indicates how incorporation of HFE in combination with OOP distortion induces and stabilizes highspin solutions on Fe center (analysis of the projected density of states for the flat and OOP pyrrolic models can be found in SI). We can surmise that even though the OOP phase is less stable and would not be thermodynamically viable compared to the flat phase, it may be stabilized by binding axial ligands or external constraints such as curving the graphene matrix in a carbon nanotube. Furthermore, our GGA+U results confirm a 2.3% decrease in work function in the 6×6 OOP compared to the flat structure, which could influence chemical reactivity between iron and adsorbates.