I have re-checked several times the following calculations and I don't know if I am missing something.
Show that the equation of motion $$ i \frac{d}{dt} \hat{\Psi}(\vec{r},t) = [\frac{1}{2m}[i\nabla + \frac{q}{c}\vec{A}(\vec{r},t)]^2 + qV(\vec{r},t)]\hat{\Psi}(\vec{r},t), \quad(1) $$ is invariant under the gauge transformation $$ \vec{A'}(\vec{r},t) = \vec{A}(\vec{r},t) + \nabla\Lambda(\vec{r},t)\\ V'(\vec{r},t) = V(\vec{r},t) -\frac{1}{c}\frac{\partial}{\partial t}\Lambda(\vec{r},t)\\ \hat{\Psi}'(\vec{r},t) = \hat{\Psi}(\vec{r},t)e^{iq\Lambda(\vec{r},t)/c}. $$
I will concentrate on the first term of equation (1) since is the one that I am stuck. In the new coordinates, the first term gives
$$ [i\nabla + \frac{q}{c}\vec{A'}(\vec{r},t)]^2 = -\nabla^2 + \frac{iq}{c}(\nabla \bullet\vec{A} + \nabla^2\Lambda) + \frac{iq}{c}(\vec{A}\bullet\nabla + \nabla\Lambda\bullet\nabla) + \frac{q^2}{c^2}(A^2+\vec{A}\bullet\nabla\Lambda + \nabla\Lambda\bullet\vec{A} + (\nabla\Lambda)^2). $$
Now, upon application of the operators to the new field operator $\hat{\Psi}'(\vec{r},t)$, I obtain
$$ -\nabla^2(\hat{\Psi} e^{iq\Lambda/c}) = [-\nabla^2\hat{\Psi}-\frac{2iq}{c}\nabla\Lambda\bullet\nabla\hat{\Psi} + \frac{q^2}{c^2}\hat{\Psi}(\nabla\Lambda)^2 - \frac{iq}{c}\hat{\Psi}\nabla^2\Lambda]e^{iq\Lambda/c}; $$
$$ \frac{iq}{c}(\nabla\bullet\vec{A} + \nabla^2\Lambda)\hat{\Psi}e^{iq\Lambda/c}; $$
$$ \frac{iq}{c}(\vec{A}\bullet\nabla + \nabla\Lambda\bullet\nabla)\hat{\Psi}e^{iq\Lambda/c} = \frac{iq}{c}(\vec{A}\bullet\nabla\hat{\Psi} + \nabla\Lambda\bullet\nabla\hat{\Psi})e^{iq\Lambda/c} -\frac{q^2}{c^2}(\vec{A}\bullet\nabla\Lambda + (\nabla\Lambda)^2)\hat{\Psi}e^{iq\Lambda/c}; $$
and finally
$$ \frac{q^2}{c^2}(A^2+\vec{A}\bullet\nabla\Lambda + \nabla\Lambda\bullet\vec{A} + (\nabla\Lambda)^2)\hat{\Psi}e^{iq\Lambda/c}. $$
Now, the trouble is that summing up this four equations, I don't get
$$ [i\nabla + \frac{q}{c}\vec{A}(\vec{r},t)]^2. $$
Any idea what went wrong?
