  Wannier functions have become a powerful tool in the electronic structure
calculations of extended systems. The generalized Pipek-Mezey Wannier functions
exhibit appealing characteristics (e.g., reaching an optimal localization and
the separation of the $\sigma$-$\pi$ orbitals) when compared with other
schemes. However, when applied to giant nanoscale systems, the orbital
localization suffers from a large computational cost overhead when one is
interested in localized states in a small fragment of the system. Herein we
present a swift, efficient, and robust approach for obtaining regionally
localized orbitals of a subsystem within the generalized Pipek-Mezey scheme.
The proposed algorithm introduces a reduced workspace and sequentially exhausts
the entire orbital space until the convergence of the localization functional.
It tackles systems with $\sim$10000 electrons within 0.5 hours with no loss in
localization quality compared to the traditional approach. Regionally localized
orbitals with a higher extent of localization are obtained via judiciously
extending the subsystem's size. Exemplifying on large bulk and a 4-nm wide slab
of diamond with NV$^-$ center, we demonstrate the methodology and discuss how
the choice of the localization region affects the excitation energy of the
defect. Furthermore, we show how the sequential algorithm is easily extended to
stochastic methodologies that do not provide individual single-particle
eigenstates. It is thus a promising tool to obtain regionally localized states
for solving the electronic structure problems of a subsystem embedded in giant
condensed systems.
