This might already be implemented, however... many applications require the evaluation of atomic overlap matrices

$$ S_{\mu \nu}^{A} = \int w_A({\bf r}) \chi_\mu ({\bf r}) \chi_\nu({\bf r}) {\rm d}^3 r $$

where $A$ is the atom and $\mu$ and $\nu$ are basis function indices.

An example use case is our generalized Pipek-Mezey orbital localization method, which replaces the original ill-defined Mulliken (or Löwdin) charges with a variety of mathematically well-defined partial charge estimates. It turns out that the localized orbitals are remarkably insensitive to the partial charge method, which can thereby be chosen by computational convenience, such as the Becke charges defined by the above overlap matrices

$$Q_{ij}^{A} = C_{\mu i} S_{\mu \nu} C_{\nu j}$$

We have also extended this method to forming generalized Pipek-Mezey Wannier functions