The orbital magnetization in periodic solids has been nicely described by the so-called modern theory of magnetization [1,2,3,4].
$$\tag{1}
\mathcal{M}_{orb} = -\frac{1}{2} \Im \sum_{n,\mathbf{k}} f_{n\mathbf{k}} w_{\mathbf{k}} \langle \partial_\mathbf{k} u_{n\mathbf{k}} | \times (\mathcal{H}_{\mathbf{k}}+\epsilon_{\mathbf{k}}-2\mu)|\partial_\mathbf{k} u_{n\mathbf{k}} \rangle.
$$
where $f_{n\mathbf{k}}$ is the Fermi occupation of band $n$ and wave vector $\mathbf{k}$. $w_{\mathbf{k}}$ is the k-point weight. $u_{n\mathbf{k}}$ is the cell periodic part of the Bloch function $\psi_{n\mathbf{k}}=e^{i\mathbf{k \cdot r}}u_{n\mathbf{k}}$ and $\mathcal{H}_{\mathbf{k}}u_{n\mathbf{k}}=\epsilon_{\mathbf{k}}u_{n\mathbf{k}}$. $\mu$ is the chemical potential. $\times$ is cross product.
I'm wondering if we could calculate $\mathcal{M}_{orb}$ in the following way (in atomic units):
\begin{align}
\mathcal{M}_{orb} &= - \frac{1}{2} \sum_{n,\mathbf{k}} f_{n\mathbf{k}} w_{\mathbf{k}} \langle \psi_{n\mathbf{k}} | \mathbf{r} \times \mathbf{v} | \psi_{n\mathbf{k}} \rangle \tag{2}\\
&= - \frac{1}{2} \sum_{n,\mathbf{k}} f_{n\mathbf{k}} w_{\mathbf{k}} \langle u_{n\mathbf{k}} | \mathbf{r} \times \mathbf{v_k} | u_{n\mathbf{k}} \rangle \tag{3}\\
&= - \frac{1}{2} \sum_{n,\mathbf{k}} f_{n\mathbf{k}} w_{\mathbf{k}} \sum_m \langle u_{n\mathbf{k}} | \mathbf{r} |u_{m\mathbf{k}}\rangle \times \langle u_{m\mathbf{k}} |\mathbf{v_k} | u_{n\mathbf{k}} \rangle \tag{4}\\
&={\small - \frac{1}{2} \sum_{n,\mathbf{k}} f_{n\mathbf{k}} w_{\mathbf{k}} \left( \mathcal{A}_{nn}(\mathbf{k}) \times \partial_\mathbf{k} \epsilon_\mathbf{k} + \frac{i}{\hbar}\sum_{m} \mathcal{A}_{nm}(\mathbf{k}) \times \mathcal{A}_{mn}(\mathbf{k})(\epsilon_{m\mathbf{k}}-\epsilon_{n\mathbf{k}}) \right)}
\tag{5}
\end{align}
where I have used in Eq. (4) the completeness relation:
$$\sum_m |u_{mk}\rangle \langle u_{mk}|=1\tag{6},$$ with $m$ running over all the states (which should be subjected to a convergence study for real calculations) and in Eq. (5) the relation: $$i\hbar \langle u_{m\mathbf{k}} |\mathbf{v_k} | u_{n\mathbf{k}} \rangle =\mathcal{A}_{mn}(\mathbf{k}) (\epsilon_{n\mathbf{k}}-\epsilon_{m\mathbf{k}})\tag{7}.$$ with $\mathcal{A}_{mn}(\mathbf{k})=i \langle u_{m\mathbf{k}} |\nabla_\mathbf{k} | u_{n\mathbf{k}} \rangle$. So the calculation of the orbital magnetization is reduced to evaluate the matrix elements of the Berry connection and band gradient.
