Gauge Invariance of Schrodinger Equation I was working out the math of Gauge Invariance of Schrodinger Equation.
Gauge transformations on my $\psi$, and the fields $\vec{A}$, and $A_0$ are
\begin{eqnarray*}
\psi(x)&\rightarrow& e^{i\alpha(x)}  \psi(x)\\
 \vec A &\rightarrow& \vec A + \frac{1}{q} \vec\nabla \alpha(x,t)\\
 A_0 &\rightarrow& A_0 - \frac{1}{q} \partial_0 \alpha(x,t)
\end{eqnarray*}

I need to show that 
  \begin{eqnarray*}
i \frac{\partial \psi'}{\partial t}&=& H' \psi'\\
\implies i \frac{\partial \psi}{\partial t}&=& H \psi
\end{eqnarray*}
  where my Hamiltonian is given by
  \begin{eqnarray*}
H=\frac{1}{2m}(p-qA)^2 + qA_0
\end{eqnarray*}

First I try the kinetic energy term, 
\begin{eqnarray*}
(p-qA')^2 \psi' &=& \left(-i \vec\nabla -qA-\vec\nabla\alpha\right)^2 (e^{i\alpha} \psi)\\
&=&\left(i \vec\nabla +qA+\vec\nabla\alpha\right)\left(i \vec\nabla +qA+\vec\nabla\alpha\right) (e^{i\alpha} \psi)\\
&=& \left(i \vec\nabla +qA+\vec\nabla\alpha\right)
\left[ ie^{i\alpha}\vec\nabla\psi
+qAe^{i\alpha}\psi
\right]\\
&=&   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- \vec\nabla
\left( e^{i\alpha}\vec\nabla\psi\right)
+i \vec\nabla\left(qAe^{i\alpha}\psi
\right)
+\left(qA+\vec\nabla\alpha\right) \left[ ie^{i\alpha}\vec\nabla\psi +qAe^{i\alpha}\psi \right]\\
\end{eqnarray*}
Now, I am worried about the second term in the last line above,
$i \vec\nabla\left(qAe^{i\alpha}\psi\right)$. I assume that should be 
\begin{eqnarray*}
i \vec\nabla\left(qAe^{i\alpha}\psi\right) =
iq\left(\vec\nabla A\right)e^{i\alpha}\psi
+iqA\left(\vec\nabla e^{i\alpha}\right)\psi
+iqAe^{i\alpha}\left(\vec\nabla\psi\right)
\end{eqnarray*}
My problem isn't solved unless $\vec \nabla A$ is zero. Now, is it zero?
What is it anyway? Divergence of $\vec A$ or gradient of $\vec A$?
If it's divergence, then that is zero only in electrostatics. If it is a gradient, then it has some other tensor component. 
 A: The problem is also solved if $\vec\nabla A$ is non-zero. Your last line equals
$$ \mathrm e^{\mathrm i\alpha} \left( -\Delta \psi + \mathrm i q\, (\vec\nabla \cdot \vec A)\, \psi + 2 \mathrm i q\, \vec A \cdot \vec\nabla \psi + (q\vec A)^2\, \psi\right) \;. $$
Note that this is the same as
$$ \mathrm e^{\mathrm i\alpha}\, ( \vec p - q\vec A ) ( \vec p - q\vec A )\, \psi \;. $$
Edit to respond to comment.
You already did everything correctly in your question, $\vec p - q\vec A$ acts as follows on a wave function $\psi(x)$:
$$ (( \vec p - q\vec A ) \psi)(x) = -\mathrm i\hbar\, \vec \nabla \psi(x) - q\vec A(x)\, \psi(x) $$
If you act again from the left with $\vec p - q\vec A$, you just take the inner product.
Thinking about this in terms of gradients, divergences etc is possible but not too helpful, it's better to just think in components:
$$ (( \vec p - q\vec A )^2 \psi)(x) = \sum_{i=1}^3 \bigl( -\mathrm i\hbar\, \partial_i - q A_i(x) \bigr) \bigl( -\mathrm i\hbar\, \partial_i \psi(x) - q A_i(x)\, \psi(x) \bigr) $$
