The Gauge invariant momentum operator is said to be $\tilde{\boldsymbol{P}}=\boldsymbol{P}-q\nabla\Lambda(\boldsymbol{R})$, where $\boldsymbol{R}$ is the position operator and $\Lambda(\boldsymbol{R})$ is a real function.
The given unitary transformation is $U=e^{iq\Lambda(\boldsymbol{R})/\hbar}$. So, in order to show the form of $\tilde{\boldsymbol{P}}$ I need to calculate:
$$\tilde{\boldsymbol{P}}=U\boldsymbol{P}U^{\dagger}=e^{iq\Lambda(\boldsymbol{R})/\hbar}\boldsymbol{P}e^{-iq\Lambda(\boldsymbol{R})/\hbar}.$$
I think I can proceed using a Taylor expansion of $U$, so:
$$\begin{align}\tilde{\boldsymbol{P}}&=\bigg(1+\frac{i}{\hbar}q\Lambda(\boldsymbol{R})+\ldots\bigg)\boldsymbol{P}\bigg(1-\frac{i}{\hbar}q\Lambda(\boldsymbol{R})+\ldots\bigg)\\ &=\boldsymbol{P}-\frac{i}{\hbar}q\boldsymbol{P}\Lambda(\boldsymbol{R})+\frac{i}{\hbar}q\Lambda(\boldsymbol{R})\boldsymbol{P}+\frac{q^2}{\hbar^2}\Lambda(\boldsymbol{R})\boldsymbol{P}\Lambda(\boldsymbol{R})+\ldots\\ &\approx \boldsymbol{P}-\frac{i}{\hbar}q\boldsymbol{P}\Lambda(\boldsymbol{R})+\frac{i}{\hbar}q\Lambda(\boldsymbol{R})\boldsymbol{P}, \end{align} $$
where, in the last approximation I have used that $q\ll1$. If I instroduce $\boldsymbol{P}=-i\hbar\nabla$ in the last step, I obtain the desired form of $\tilde{\boldsymbol{P}}$ from the first two terms, but the last term is additional and is equivalent to $\Lambda(\boldsymbol{R})\nabla$.
Why am I obtaining this additional term? Maybe I'm just not following the right procedure. Any Suggestions?