  We present a first-principles scheme that allows the orbital magnetization of
a magnetic crystal to be evaluated accurately and efficiently even in the
presence of complex Fermi surfaces. Starting from an initial
electronic-structure calculation with a coarse ab initio k-point mesh,
maximally localized Wannier functions are constructed and used to interpolate
the necessary k-space quantities on a fine mesh, in parallel to a
previously-developed formalism for the anomalous Hall conductivity [X.Wang, J.
Yates, I. Souza, and D. Vanderbilt, Phys. Rev. B 74, 195118 (2006)]. We
formulate our new approach in a manifestly gauge-invariant manner, expressing
the orbital magnetization in terms of traces over matrices in Wannier space.
Since only a few (e.g., of the order of 20) Wannier functions are typically
needed to describe the occupied and partially occupied bands, these Wannier
matrices are small, which makes the interpolation itself very efficient. The
method has been used to calculate the orbital magnetization of bcc Fe, hcp Co,
and fcc Ni. Unlike an approximate calculation based on integrating orbital
currents inside atomic spheres, our results nicely reproduce the experimentally
measured ordering of the orbital magnetization in these three materials.
