Computational method: All the calculations were performed by applying the spin-polarized density functional theory (DFT) method implemented in the Vienna Ab initio Simulation Package (VASP). The cut-off energy was set to be 500 eV for the 2×2×2 KFeFe-PBA supercells. The DFT-D3 method was used to describe Van der Waals (vdW) interactions. The DFT + U method was used to correctly characterize the localization of the Fe d-electrons in the KFeFe-PBA supercells. To differentiate high-spin (HS) and low-spin (LS) Fe ions in the structures, the values of U-J for LS-Fe and HS-Fe were set independently to be 3.0 and 7.0 eV in this work, respectively. During the geometry optimization, the convergence criteria of the force and energy change for all atomic structures were set to be 0.05 eV Å -1 and 10 -5 eV, respectively. The minimum energy pathway search of ion migration was conducted with the climbing image nudged elastic band (CI-NEB) method. The intercalation energy (𝐸 𝑖 ) of the cations was given by