Derive probability current density - factors of 2 discrepancy To derive the probability current density for a particle
in an electromagnetic field, we calculate
$\dfrac{\partial \rho}{\partial t}
=
\dfrac{\partial}{\partial t} (\Psi^* \Psi)
=
\dfrac{\partial \Psi^*}{\partial t} \Psi
+
\Psi^* \dfrac{\partial \Psi}{\partial t}$
$H$ is, if we substitute in $-i\hbar \nabla$ for $\vec{p}$:
$H = \frac{1}{2m}(\vec{p} - \frac{q}{c} \mathbf{A}) \cdot
                 (\vec{p} - \frac{q}{c} \mathbf{A})   + q\phi
   = \frac{1}{2m}(-i\hbar \nabla - \frac{q}{c} \mathbf{A}) \cdot
                 (-i\hbar \nabla - \frac{q}{c} \mathbf{A})   + q\phi
   = \frac{1}{2m}(i\hbar \nabla + \frac{q}{c} \mathbf{A}) \cdot
                 (i\hbar \nabla + \frac{q}{c} \mathbf{A})   + q\phi$
Schrödinger equation and its complex conjugate:
$\dfrac {\partial \Psi}{\partial t} = \dfrac{H\Psi}{i\hbar}$
$\dfrac {\partial \Psi^*}{\partial t} = \dfrac{(H \Psi)^*}{-i\hbar}$
Substitute in:
$\dfrac {\partial \rho}{\partial t}
=
\dfrac{-1}{i\hbar} [(H\Psi)^* \Psi - \Psi^* (H\Psi)]$
Substitute in $H$:
$\dfrac {\partial \rho}{\partial t}
= \dfrac{-1}{i\hbar}
 \{ ( [\frac{1}{2m}(+i\hbar \nabla + \frac{q}{c} \mathbf{A})
        \cdot      (+i\hbar \nabla + \frac{q}{c} \mathbf{A}) + q\phi]\Psi )^* \Psi \\
 - \Psi^*( [\frac{1}{2m}(+i\hbar \nabla + \frac{q}{c} \mathbf{A})
             \cdot      (+i\hbar \nabla + \frac{q}{c} \mathbf{A}) + q\phi] \Psi)
 \}$
Apply complex conjugate:
$\dfrac {\partial \rho}{\partial t}
=\dfrac{-1}{i\hbar}
 \{ ( [\frac{1}{2m}(-i\hbar \nabla + \frac{q}{c} \mathbf{A})
        \cdot      (-i\hbar \nabla + \frac{q}{c} \mathbf{A}) + q\phi]\Psi^* ) \Psi \\
 - \Psi^*( [\frac{1}{2m}(+i\hbar \nabla + \frac{q}{c} \mathbf{A})
             \cdot      (+i\hbar \nabla + \frac{q}{c} \mathbf{A}) + q\phi] \Psi)
 \}$
FOIL:
$\dfrac {\partial \rho}{\partial t}=\dfrac{-1}{i\hbar}
 \{ ( [\frac{1}{2m}(i\hbar i\hbar \nabla^2 + (-i\hbar) \nabla \cdot (\frac{q}{c} \mathbf{A})
    + (\frac{q}{c} \mathbf{A}) \cdot (-i\hbar) \nabla + \frac{q^2}{c^2} \mathbf{A}^2)
    + q\phi]\Psi^* ) \Psi \\
 - \Psi^*( [\frac{1}{2m}(i\hbar i\hbar \nabla^2 + i\hbar \nabla \cdot (\frac{q}{c} \mathbf{A})
    + (\frac{q}{c} \mathbf{A}) \cdot (i\hbar \nabla) + \frac{q^2}{c^2} \mathbf{A}^2)
    + q\phi] \Psi)
 \}$
Multiply everything out:
$\dfrac {\partial \rho}{\partial t}
=\frac{-i\hbar}{2m}(\nabla^2 \Psi^*) \Psi
+ \frac{1}{2m} (\nabla \cdot \frac{q}{c} \mathbf{A}) \Psi^* \Psi
+ \frac{1}{2m} (\frac{q}{c} \mathbf{A}) \cdot (\nabla \Psi^*) \Psi
+ \frac{-1}{i\hbar} \frac{1}{2m} \frac{q^2}{c^2} \mathbf{A}^2 \Psi^* \Psi \\
+ \frac{-1}{i\hbar} \frac{1}{2m} q \phi \Psi^* \Psi \\
+ (\Psi^*) \frac{1}{2m} (i\hbar)(\nabla^2 \Psi)
+ (\Psi^*) \frac{1}{2m} \nabla \cdot (\frac{q}{c} \mathbf{A}) \Psi
+ (\Psi^*) \frac{1}{2m} (\frac{q}{c} \mathbf{A}) \cdot (\nabla \Psi)
+ \frac{1}{i\hbar} (\Psi^*) \frac{1}{2m} \frac{q^2}{c^2} \mathbf{A}^2 \Psi \\
+ \frac{1}{i\hbar} \frac{1}{2m}(\Psi^*)q\phi \Psi
$
The terms containing  $\phi$  and  $\frac{q^2}{c^2} \mathbf{A}^2$  cancel
and there's a fact that   $\Psi \nabla^2 \Psi^* - \Psi^* \nabla^2 \Psi
= \nabla \cdot(\Psi \nabla \Psi^* - \Psi^* \nabla \Psi)$, so
$\dfrac {\partial \rho}{\partial t}
= \frac{-i\hbar}{2m} \nabla \cdot (\Psi \nabla \Psi^* - \Psi^* \nabla \Psi) \\
+ \frac{1}{2m} (\nabla \cdot \frac{q}{c} \mathbf{A}) \Psi^* \Psi
+ \frac{1}{2m} (\frac{q}{c} \mathbf{A}) \cdot (\nabla \Psi^*) \Psi
+ (\Psi^*) \frac{1}{2m} \nabla \cdot (\frac{q}{c} \mathbf{A}) \Psi
+ (\Psi^*) \frac{1}{2m}  (\frac{q}{c} \mathbf{A}) \cdot (\nabla \Psi)$
Note that of the 5 terms, the 2nd and 4th are the same, so
(1) $\dfrac {\partial \rho}{\partial t}
= \frac{-i\hbar}{2m} \nabla \cdot (\Psi \nabla \Psi^* - \Psi^* \nabla \Psi) \\
+ \dfrac{q}{mc} (\nabla \cdot \mathbf{A}) \Psi^* \Psi
+ \dfrac{q}{2mc} \mathbf{A} \cdot (\Psi \nabla \Psi^*)
+ \dfrac{q}{2mc} \mathbf{A} \cdot (\Psi^* \nabla \Psi)$
https://en.wikipedia.org/wiki/Probability_current tells us
that the final result should be $\dfrac{\partial \rho}{\partial t}
 = - \nabla \cdot \mathbf{j}$ and that
$\mathbf{j} = \dfrac{1}{2m} [(\Psi^* \mathbf{\hat{p}} \Psi
- \Psi \mathbf{\hat{p}} \Psi^* ) - 2 \frac{q}{c} \mathbf{A} |\Psi|^2]$
or using $\mathbf{\hat{p}} = -i\hbar \nabla$,
$\mathbf{j} = \dfrac{1}{2m} [(\Psi^* (-i\hbar \nabla) \Psi
- \Psi (-i\hbar \nabla) \Psi^* ) - 2 \frac{q}{c} \mathbf{A} |\Psi|^2]$
$\mathbf{j} = \dfrac{-i\hbar}{2m} (\Psi^* \nabla \Psi
- \Psi \nabla \Psi^* ) - \dfrac{1}{2m} 2 \frac{q}{c} \mathbf{A} |\Psi|^2$
$\mathbf{j} = \dfrac{-i\hbar}{2m} (\Psi^* \nabla \Psi
- \Psi \nabla \Psi^* ) - \dfrac{q}{mc} \mathbf{A} |\Psi|^2$
Applying $\dfrac{\partial \rho}{\partial t} = - \nabla \cdot \mathbf{j}$,
$\dfrac{\partial \rho}{\partial t}
= - \nabla \cdot [\dfrac{-i\hbar}{2m} (\Psi^* \nabla \Psi
- \Psi \nabla \Psi^* ) - \dfrac{q}{mc} \mathbf{A} |\Psi|^2]$
$\dfrac{\partial \rho}{\partial t}
= \dfrac{i\hbar}{2m} \nabla \cdot (\Psi^* \nabla \Psi
- \Psi \nabla \Psi^* ) + \dfrac{q}{mc} \nabla \cdot
(\mathbf{A} |\Psi|^2)$
Apply an identity about $\nabla$ operating on a scalar times a vector
at https://en.wikipedia.org/wiki/Vector_calculus_identities
$\nabla \cdot (\phi \mathbf{B}) =
\mathbf{B} \cdot \nabla \phi + 
\phi (\nabla \cdot \mathbf{B})$   (I changed the letters),
$\dfrac{\partial \rho}{\partial t}
= \dfrac{i\hbar}{2m} \nabla \cdot (\Psi^* \nabla \Psi
- \Psi \nabla \Psi^* )
+ \dfrac{q}{mc} (\mathbf{A} \cdot \nabla (\Psi^* \Psi)
+ (\Psi^* \Psi) (\nabla \cdot \mathbf{A}))$
Product rule:
$\dfrac{\partial \rho}{\partial t}
= \dfrac{i\hbar}{2m} \nabla \cdot (\Psi^* \nabla \Psi
- \Psi \nabla \Psi^* )
+ \dfrac{q}{mc} (\mathbf{A} \cdot (\Psi^* \nabla \Psi)
               + \mathbf{A} \cdot (\Psi \nabla \Psi^*)
+ (\Psi^* \Psi) (\nabla \cdot \mathbf{A}))$
$\dfrac{\partial \rho}{\partial t}
= \dfrac{i\hbar}{2m} \nabla \cdot (\Psi^* \nabla \Psi
- \Psi \nabla \Psi^* )
+ \dfrac{q}{mc} \mathbf{A} \cdot (\Psi^* \nabla \Psi)
+ \dfrac{q}{mc} \mathbf{A} \cdot (\Psi \nabla \Psi^*)
+ \dfrac{q}{mc} (\Psi^* \Psi) (\nabla \cdot \mathbf{A})$
Rearrange some terms so that we can compare with equation (1):
$\dfrac{\partial \rho}{\partial t}
= \dfrac{i\hbar}{2m} \nabla \cdot (\Psi^* \nabla \Psi
- \Psi \nabla \Psi^* )
+ \dfrac{q}{mc} (\Psi^* \Psi) (\nabla \cdot \mathbf{A})
+ \dfrac{q}{mc} \mathbf{A} \cdot (\Psi \nabla \Psi^*)
+ \dfrac{q}{mc} \mathbf{A} \cdot (\Psi^* \nabla \Psi)$
Now we're very close to equation (1),
(1) $\dfrac {\partial \rho}{\partial t}
= \frac{-i\hbar}{2m} \nabla \cdot (\Psi \nabla \Psi^* - \Psi^* \nabla \Psi) \\
+ \dfrac{q}{mc} (\nabla \cdot \mathbf{A}) \Psi^* \Psi
+ \dfrac{q}{2mc} \mathbf{A} \cdot (\Psi \nabla \Psi^*)
+ \dfrac{q}{2mc} \mathbf{A} \cdot (\Psi^* \nabla \Psi)$
but are off by some factors of $2$.
Do you see an error in my steps?
 A: 
Do you see an error in my steps?

Yes, I see a conceptual error. So I'll talk about that.

$\dfrac{\partial \rho}{\partial t}
=
\dfrac{\partial}{\partial t} (\Psi^* \Psi)
=
\dfrac{\partial \Psi^*}{\partial t} \Psi
+
\Psi^* \dfrac{\partial \Psi}{\partial t}$

should be
$\dfrac{\partial \rho}{\partial t}
=
\dfrac{\partial}{\partial t} (\Psi^* \Psi)
=
\Psi \dfrac{\partial }{\partial t} 
\Psi^*
+
\Psi^* \dfrac{\partial }{\partial t}\Psi$
because $\partial_t$ is an operator and so we should only have things to the right of it that it acts on. Because you might remember that now but later you will forget and you will then make a mistake. Plus there are going to just be too many parenthesis anyway.
Now look at

$H = \frac{1}{2m}(\vec{p} - \frac{q}{c} \mathbf{A}) \cdot
                 (\vec{p} - \frac{q}{c} \mathbf{A})   + q\phi
   = \frac{1}{2m}(-i\hbar \nabla - \frac{q}{c} \mathbf{A}) \cdot
                 (-i\hbar \nabla - \frac{q}{c} \mathbf{A})   + q\phi
$ 

which equals $$\frac{1}{2m}\sum_j\left((i\hbar \partial_j +\frac{q}{c} A_j)(i\hbar \partial_j +\frac{q}{c} A_j)\right)  + q\phi.$$
And these are operators the one on the right acts on things then the next acts on the outcome of that and so on. Think of them like matrices. Don't imagine any vectors or dot products don't start thinking something (anything) is a scalar and don't start changing the order of anything. If you couldn't justify it if it were a matrix, then don't do it. Just treat them all like matrices. Even a scalar field like $A_j$ first it has to multiply something then a derivative has to act on that product (with the product rule) otherwise you are doing it wrong.
