I’m trying to derive the equation from the Hamiltonian
$$H=\frac{1}{m}[\vec{\sigma} \cdot (\vec p-q\vec A)]^2+ q \phi$$
and
$$[\vec{\sigma} \cdot (\vec p-q\vec A)]^2= (\vec p-q\vec A)^2-iq \vec{\sigma} \cdot \nabla \times (\vec p-q\vec A).$$
I don’t know how the term $\nabla \times \vec p$ vanish in the result...
Using relation $\sigma_i\sigma_j=\delta_{ij}+i\epsilon_{ijk}\sigma_k$ and denotation $P_i=p_i-eA_i$, we can write $$(\sigma_iP_i)(\sigma_jP_j)=\delta_{ij}P_iP_j+i\epsilon_{ijk}P_iP_j\sigma_k,$$ where the last term explicitly is $$P_iP_j=p_ip_j-ep_iA_j-eA_ip_j+e^2A_iA_j.$$ Due to $\epsilon_{ijk}p_ip_j=0$ and $\epsilon_{ijk}A_iA_j=0$, you should find $$-i\epsilon_{ijk}p_iA_j\sigma_k-i\epsilon_{ijk}A_ip_j\sigma_k\rightarrow+(\nabla\times {\bf A})\cdot\sigma=\sigma\cdot{\bf H}.$$
In your notation, you should see that $$\nabla\times{\bf p}=\epsilon_{ijk}\partial_jp_k=-i\epsilon_{ijk}\partial_j\partial_k,$$ which goes to zero due to the fact that ${\bf v}\times{\bf v}=0$, where ${\bf v}$ is arbitrary vector.
