Attempt at proving $-i langle{u_n|nabla_k u_n}rangle=-dfrac{i}{2}tr[v^dagger(k)nabla_k v(k) ]$ from Kane and Fu's paper

I am trying to prove result (3.4) of the following paper:
http://li.mit.edu/S/2d/Paper/Fu07Kane.pdf
namely, that $$-i langle{u_n|nabla_k u_n}rangle=-dfrac{i}{2}tr[v^dagger(k)nabla_k v(k) ]$$
where I am using the Einstein summation convention and $v$ is given by
$$v_{nm}(k)=langle{u_n|PTheta|u_mrangle}$$
where $Theta$ is the time reversal operator and $P$ is the space inversion operator, and all the $langle{u_n|}$ are evaluated at the same $k$ (this is a vector).

My attempt is the following:
I first write $-ilangle{u_n|nabla_k u_nrangle}=-langle{u_n|r|u_nrangle}$ because I know that both $Theta$ takes $rmapsto r$ and $P$ takes $rmapsto -r$. I will also use that $v v^+=v^dagger v=1$, $P^dagger P=1$ and $langle{Theta psi |Theta phirangle}=langle{phi|psirangle}=langlepsi|phirangle^*$.
So, here it goes:
$$-langle{ u_n |r|u_n rangle}=-langle{ PTheta r u_n | PTheta u_n rangle}=langle{ r PTheta u_n | PTheta u_n rangle}=langle{ r PTheta u_n | u_m rangle}langle{u_m |PTheta u_n rangle}=v_{mn}langle{ r PTheta u_n | u_m rangle}$$
I now use the chain rule (recalling that the $k-$space representation of $r$ is a derivative with respect to $k$) to write for example $rlangle{psi|phirangle}=langle{r psi|phirangle}+langle{psi|r phirangle}$. Doing so, gives:
$$langle{ u_n |r|u_n rangle}=v_{mn}left(r langle{PTheta u_n|u_mrangle}-langle{PTheta u_n|r u_mrangle} right)$$
I now insert the resolution of the identity in the last term:
$$-langle{ u_n |r|u_n rangle}=v_{mn}left( r v^*_{mn}-langle{PTheta u_n|u_s rangle}langle{u_s|r u_mrangle} right)=v_{mn}r v^*_{mn}-underbrace{v_{mn}v^*_{sn}}_{delta_{ms}}langle{u_s|r u_mrangle}$$
$$=tr[v r v^dagger]-langle{u_m|r|u_mrangle}=tr[v r v^dagger]-langle{u_n|r|u_nrangle}$$
Which is clearly wrong!

I suspect that a mistake is that $$-langle{ PTheta r u_n | PTheta u_n rangle}=langle{ r PTheta u_n | PTheta u_n rangle}$$, but even if I don’t change the sign here, the last result picks up a minus sign, which gives that $$-i langle{u_n|nabla_k u_n}rangle=-dfrac{i}{2}tr[v(k)nabla_k v^dagger(k) ]=dfrac{i}{2}tr[v^dagger(k)nabla_k v(k) ]$$
which differs by a sign from the result I should have gotten; so, even if that sign change when I (anti)commute $r$ past $PTheta$ doesn’t fix things.
Any help is appreciated, thanks!