## 2 Methodology  
The calculations were performed using the Vienna Ab Initio Simulation Package (VASP) developed at the Universitat Wien [23; 24]. Our previous research showed that the strongly constrained and appropriately normed (SCAN) functional [25] correctly reproduces the electrical and magnetic properties of the YSr\({}_{2}\)Cu\({}_{2}\)FeO\({}_{7}\) and YSr\({}_{2}\)Cu\({}_{2}\)FeO\({}_{8}\) end members. Therefore, it was employed to study the intermediate compositions of the YSr\({}_{2}\)Cu\({}_{2}\)FeO\({}_{7+\delta}\) system. The interaction of core electrons with the nuclei is described by the projector augmented wave (PAW) method [26], with 2\(s\)2\(p\) of O, 3\(d\)4\(s\) of Fe, 3\(d\)4\(p\)5s of Sr, and 4\(s\)4\(p\)5s4\(d\) of Y treated as valence electrons. The energy cutoff for the plane wave basis set was fixed at a constant value of 600 eV throughout the calculations. The integration in the Brillouin zone is performed on appropriate sets of k-points, determined by the Monkhorst-Pack scheme. The k-point meshes were set at 10 \(\times\) 10 \(\times\) 4 using a Gaussian smearing parameter of 0.05 eV. For the density of states (DOS) calculations, the tetrahedron method with Bloch corrections was used [27]. Self-consistency was achieved with a tolerance in total energy of 1 \(\times\) 10\({}^{-4}\) eV for geometry optimization and 1 \(\times\) 10\({}^{-6}\) eV for DOS calculations. To replicate the oxygen non-stoichiometry, the end-term YSr\({}_{2}\)Cu\({}_{2}\)FeO\({}_{8}\) (S.G. P4/_mmm_) is considered as the non-defective phase. Following our previous work, the starting positions for YSr\({}_{2}\)Cu\({}_{2}\)FeO\({}_{8}\) are taken from (ICSD file 11514, see Figure 1), considering a 2 \(\times\) 2 \(\times\) 1 superstructure [20]. Hence, the lattice parameters of the idealized model of YSr\({}_{2}\)Cu\({}_{2}\)FeO\({}_{8}\) are 2a\({}_{\text{p}}\)\(\times\) 2a\({}_{\text{p}}\)\(\times\) 3a\({}_{\text{p}}\). To model the intermediate oxygen contents (\(\delta\) = 0.2, 0.5, and 0.75), oxygen vacancies are accommodated in the O3 position of the FeO\({}_{1+\delta}\) layers. The construction of oxygen/vacancies configurations was carried out using the Pymatgen Python library (v.2023.9.10) [28]. All the models are considered in a ferromagnetic configuration. The Ewald summation method [29] was employed to choose, out of those generated, the most stable configurations for DFT calculations. The local magnetic moments and charges of each atom were obtained from the Bader charge analysis [30], which was performed on the charge density files [31] using the Pymatgen Python library [28]. For the ab initio molecular dynamics (AIMD) calculations, the generalized gradient approximation exchange and correlation functional developed by [PERSON] (PBE) [32] was used. Sampling of the electronic structure was made at the gamma point only. Calculations at 1000, 1200, 1400, 1600, and 2000 K in the NVT ensemble were performed for YSr\({}_{2}\)Cu\({}_{2}\)FeO\({}_{7.5}\). First, the temperature was increased during 2000 steps, and it was then held constant for 12,000 successive steps, employing the Nose-Hover thermostat [33]. A cutoff energy of 400 eV was used, along with a time step of 2 fs for all calculations. The data analysis was performed using the Pymatgen-diffusion [34] add-on for the Pymatgen [28] package. Ionic transport was monitored through mean square displacement (MSD). Studying the MSD variation of ions as a function of time, the diffusivity D can be extracted from the slopes for each simulated temperature, according to:  
\[\langle|\text{r}_{\text{i}}(\text{t})-\text{r}_{\text{i}}(0)|\rangle=6\text{D} \text{t}+\text{B} \tag{1}\]  
The activation energy (E\({}_{\text{a}}\)) of the process can be calculated from the linear fitting of the obtained data to the Arrhenius plot:  
\[\text{D}\ =\text{D}_{0}\ \text{exp}\bigg{(}-\frac{\text{E}_{\text{a}}}{\text{kT}} \bigg{)} \tag{2}\]where \(\rm D_{0}\) is the oxygen diffusivity at \(\rm T\rightarrow\infty\), that is to say, the maximum oxygen diffusivity.