  Maximally localized Wannier functions (MLWFs) are widely used to construct
first-principles tight-binding models that accurately reproduce the electronic
structure of materials. Recently, robust and automated approaches to generate
these MLWFs have emerged, leading to natural sets of atomic-like orbitals that
describe both the occupied states and the lowest-lying unoccupied ones (when
the latter can be meaningfully described by bonding/anti-bonding combinations
of localized orbitals). For many applications, it is important to instead have
MLWFs that describe only certain target manifolds separated in energy between
them -- the occupied states, the empty states, or certain groups of bands.
Here, we start from the full set of MLWFs describing simultaneously all the
target manifolds, and then mix them using a combination of parallel transport
and maximal localization to construct orthogonal sets of MLWFs that fully and
only span the desired target submanifolds. The algorithm is simple and robust,
and it is applied to some paradigmatic but non-trivial cases (the valence and
conduction bands of silicon, the top valence band of MoS$_2$, the $3d$ and
$t_{2g}$/$e_g$ bands of SrVO$_3$) and to a mid-throughput study of 77
insulators.
