Confused with the velocity operator The velocity operator is defined as the commutator of the position operator and the Hamiltonian
$$
\mathbf{v} = -\frac{i}{\hbar}[\mathbf{r},H]
$$
Say $H$ is a crystal Hamiltonian with eigenstates of Bloch functions $|n\mathbf{k}\rangle$, $H|n\mathbf{k}\rangle=E_{n\mathbf{k}}|n\mathbf{k}\rangle$,$n$ the band index and $\mathbf{k}$ the wave vector. We know, in general, the velocity of Bloch electron $\langle n\mathbf{k}|\mathbf{v}|n\mathbf{k}\rangle=1/\hbar\partial_\mathbf{k}E_{n\mathbf{k}}$ is not zero. However, if we replace $\mathbf{v}$ with the commutator
$$
\langle n\mathbf{k}|\mathbf{v}|n\mathbf{k}\rangle = -\frac{i}{\hbar}\langle n\mathbf{k}|\mathbf{r}H-H\mathbf{r}|n\mathbf{k}\rangle = -\frac{i}{\hbar}(E_{n\mathbf{k}}-E_{n\mathbf{k}})\langle n\mathbf{k}|\mathbf{r}|n\mathbf{k}\rangle = 0
$$
Am I miss something?
 A: I think that this is a problem with using the commutation relations directly on eigenstates, especially when dealing with infinite matrices. For example $[x, p] = i$ but $\langle p | [x, p] | p \rangle = (p-p)\langle p | x | p \rangle = 0$. Which all comes down to that for finite-size matrices ${\rm Tr}[A , B ] = 0$.
The "correct" way to do that is by taking a limit
$$ \lim_{p' \to p} \langle p | [x, p] | p'\rangle = \lim_{p'\to p} (p'-p)\langle p | x | p'\rangle = \lim_{p'\to p} i(p'-p)\partial_p \langle p' | p \rangle = \lim_{p'\to p}i(p'-p)\frac{\delta(p'-p)}{(p'-p)} = i\langle p | p \rangle$$
where I used $\partial_x \delta(x) = -\delta(x)/x$ but any representation of the delta function or of the $|p\rangle$ wave functions would have sufficed.
Similarly
$$ -i\langle n {\bf k} | [{\bf r}, H] | n {\bf k}'\rangle = (E_{\bf k}'-E_{\bf k})\partial_{\bf k}\langle n {\bf k} |  n {\bf k}'\rangle$$
which will result in the desired expression upon taking the limit ${\bf k}'\to {\bf k}$. You can also use the Hellman-Feynman theorem to take the limit as derivative with respect to ${\bf k}$.
