  Orbital magnetic susceptibility involves rich physics such as interband
effects despite of its conceptual simplicity. In order to appreciate the rich
physics related to the orbital magnetic susceptibility, it is essential to
derive a formula to decompose the susceptibility into the contributions from
each band. Here, we propose a scheme to perform this decomposition using the
modified Wannier functions. The derived formula nicely decomposes the
susceptibility into intraband and interband contributions, and from the other
aspect, into itinerant and local contributions. The validity of the formula is
tested in a couple of simple models. Interestingly, it is revealed that the
quality of the decomposition depends on the degree of localization of the used
Wannier functions. The formula here complements another formula using Bloch
functions, or the formula derived in the semiclassical theory, which deepens
our understanding of the orbital magnetic susceptibility and may serve as a
foundation of a better computational method. The relationship to the Berry
curvature in the present scheme is also clarified.
