The DFT calculations were performed, using version 6.2.1 of the plane-wave method based code Vienna Ab-initio Simulation Package (VASP). The projector augmented wave (PAW) approach was utilized to replace the core electrons. To address the known issue of overlocalization, the PBE functional was employed together with Dudarev's version of the Hubbard U -correction, which was applied to the iron d-states. The U values were adapted from previous studies, which indicated that PBAs are best represented with distinct U -values for each iron site. In this case, a Hubbard correction of 7 eV was assigned to the d-states of the high-spin iron within the FeN 6 octahedra, while a value of 3 eV was used for the d-states of the low-spin iron with the FeC 6 octahedra. Atomic structures were optimized until the energy difference in the electronic self-consistent field (SCF) fell lower than 10 -6 eV, and a force convergence criteria of 10 -2 eV/ Å was employed. Furthermore, all calculations used a Gamma-centered k-point grid of 5×5×5 for cubic structures and 5×5×3 for rhombohedral ones. The plane-wave cutoff was set to 550 eV and Gaussian smearing with a width of 0.1 eV was used. Additionally, the total excess majority spin electrons were held constant to ensure convergence to the correct high-spin-low-spin spin state on the iron atoms.