The convenient way is implemented in wannier90's "berry" module, where the orbital magnetization might be calculated in formalism of Berry phases. As far as I know, Wien2K has the interface with Wannier90-2.0.0, which has this module. VASP also has the interface, but with only 1.2 version, where spinor representation of Wannier functions and, consequently, "berry" module are not availiable.