Modifying the Hamiltonian when there is a presence of the Coulomb interaction Referring to the Hamiltonian of a system of free electrons,
$$
 H_0= \sum_{\sigma} \int d^3rd^3r' \psi_{\sigma}^{\dagger}(\mathbf{r})\left(- \frac{\hbar^2}{2m}\nabla^2\right)\delta(\mathbf{r}-\mathbf{r'})\psi_\sigma (\mathbf{r'})
$$
When the Coulomb interaction is turn on, we can modify this Hamiltonian by imposing
$$
\partial^\mu\rightarrow\partial^\mu + i \frac{q}{\hbar c} A^\mu
$$
I expected the new Hamiltonian to contain the term like
$$
H_{\mathrm{Coulomb}} = \frac{1}{2}\sum_{\sigma,\sigma'}\int d^3rd^3r' \psi_\sigma^\dagger(\mathbf{r})\psi_\sigma(\mathbf{r})\left( \frac{q^2}{4\pi\epsilon_0|\mathbf{r}-\mathbf{r'}|}\right)\psi_{\sigma'}^\dagger(\mathbf{r'})\psi_{\sigma'}(\mathbf{r'})
$$
However, when I substituted this to the free Hamiltonian, I could not see anyway to obtain this result at all.
 A: There are a couple of issues with the phrasing of the question:
First: Coulomb interaction only affects the temporal component of the derivative operator (i.e. you need to focus on only $\mu = 0$). The general form you wrote is useful if you consider the full electromagnetic field (i.e. in the presence of both $A_0$ and $\mathbf A$).
Second: The Peierls substitution for $\partial^0$ leads to a minimal coupling between the electric potential ($A_0$) and electron density, which is unavailable in the hamiltonian formalism. To see this the action is required. 
Within the action formalism the $A_0$ field can be integrated out, which will lead to the effective long-range interaction you have quoted. To demonstrate, I use the imaginary time ($\tau$) formalism and slightly modify the parameters in your equations to write the action,
$$S = \sum_\sigma \int d\tau  d\mathbf{r}  ~\Psi_\sigma^\dagger(\tau, \mathbf r) \left[\partial_\tau - i g A_0(\tau, \mathbf r) \right] \Psi_\sigma(\tau, \mathbf r) - \int d\tau H_0 + \frac{1}{2}\int d\tau  d\mathbf{r} |\mathbf \nabla A_0(\tau, \mathbf r)|^2.$$
Now from within the partition function $$ Z = \int DA_0 ~ D\Psi ~ D\Psi^{\dagger} ~ e^{-S},$$ one may integrate out $A_0$ (easier to do that if you Fourier transform and write $S$ in momentum space) using properties of Gaussian integration. This will result in $$S =  \sum_\sigma \int d\tau  d\mathbf{r}  ~\Psi_\sigma^\dagger(\tau, \mathbf r) \partial_\tau \Psi_\sigma(\tau, \mathbf r) - \int d\tau (H_{Coulomb} + H_0),$$ where $H_{Coulomb}$ is the equation noted in the question. Therefore, the total hamiltonian is $H = H_{Coulomb} + H_0$.
