The doped 59 defect simulation package was used for all steps in the defect modelling workflow, including structure and input file generation, calculation parsing, analysis and plotting. The standard supercell approach for computing defect formation energies was employed, and self-consistent defect/carrier concentrations were computed using the total charge neutrality condition (Section S4). For defect calculations, an 81-atom supercell was used, produced from a 3 ⇥ 3 ⇥ 3 expansion of the primitive P3 1 21 unit cell, along with a 2 ⇥ 2 ⇥ 2 k-point mesh (rounding up from the converged k-point density of 0.42 Å 1 ). The ShakeNBreak 63-65 defect structure-searching method was employed in all cases, finding many low-energy defect geometries which are missed by standard gradient relaxations of unperturbed or 'rattled' defect supercells (with DE up to 1 eV), particularly for extrinsic impurities. Further details are provided in Section S5.5.