Proof of the Gauge-invariant momentum operator 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?
 A: You should remember that these operators act on wavefunctions:
   $$UPU^+\psi=P\psi-{i\over\hbar}P\big(q\Lambda\psi\big)
   +{i\over\hbar}q\Lambda P\psi$$
to first order in $q\Lambda$. Note that in the second term in the r.h.s., $P$ acts on $\Lambda\psi$ and not only on $\Lambda$. Using the fact that
   $$P\big(\Lambda\psi\big)=(P\Lambda)\psi+\Lambda\psi$$
The last term of the first equation is therefore cancelled by the last term of the second relation. It follows that
   $$UPU^+\psi=P\psi-{iq\over\hbar}(P\Lambda)\psi
   =P\psi+q(\nabla\Lambda)\psi$$
A: If $\Lambda$ is a function of the position operator, you can't just treat $\Lambda$ as a scalar function. In fact, what you are obtaining is $\tilde{\mathbf{P}} = \mathbf{P} - i q [\mathbf{P},\Lambda(\mathbf{R})]$.
For analytic functions you can always expand $\Lambda$ as a power series in $\mathbf{R}$ and compute the commutator with $\mathbf{P}$ for each power of $\mathbf{R}$ using the canonical commutation relations and the product rule for operators. You will obtain $[\mathbf{P},\Lambda(\mathbf{R})] = - i \Lambda'(\mathbf{R})$. Using this you get the result you were looking for.
