  We introduce a new type of Wannier functions (WFs) obtained by minimizing the
conventional spread functional with a penalty term proportional to the variance
of the spread distribution. This modified Wannierisation scheme is less prone
to produce ineffective solutions featuring one or several poorly localized
orbitals, making it well suited for complex systems or high-throughput
applications. Furthermore, we propose an automatable protocol for selecting the
initial guess and determine the optimal number of bands (or equivalently WFs)
for the localization algorithm. The improved performance and robustness of the
approach is demonstrated for a diverse set of test systems including the NV
center in diamond, metal slabs with atomic adsorbates, spontaneous polarization
of ferroelectrics and 30 inorganic monolayer materials comprising both metals
and semiconductors. The methods are implemented in Python as part of the Atomic
Simulation Environment (ASE).
