We present a study of the construction and spatial properties of localized Wannier orbitals in large supercells of insulating solids using plane waves as the underlying basis. The Pipek-Mezey (PM) functional in combination with intrinsic atomic orbitals (IAOs) as projectors is employed, resulting in so-called intrinsic bond orbitals (IBOs). Independent of the bonding type and band gap, a correlation between orbital spreads and geometric properties is observed. As a result, comparable sparsity patterns of the Hartree–Fock exchange matrix are found across all considered bulk 3D materials, exhibiting covalent bonds, polar covalent bonds, and ionic bonds. Recognizing the considerable computational effort required to construct localized Wannier orbitals for large periodic simulation cells, we address the performance and scaling of different solvers for the localization problem. This includes the Broyden–Fletcher–Goldfarb–Shanno (BFGS), Conjugate-Gradient (CG), Steepest Ascent (SA), as well as the Direct Inver-sion in the Iterative Subspace (DIIS) method. Each algorithm performs a Riemannian optimization under a unitary matrix constraint, efficiently reaching the optimum in the “curved parameter space” on geodesics. The solvers have been implemented both within the VASP and as a standalone open-source software package. Furthermore, we observe that the construction of Wannier orbitals for supercells of metal oxides presents a significant challenge, requiring approximately 1 order of magnitude more iteration steps than other systems studied.