How to derive the particle current operators for the non-interacting and interacting Hubbard model ?
Hubbard Hamiltonian is given here with the interaction term: http://en.wikipedia.org/wiki/Hubbard_model
Here's an introduction : http://quest.ucdavis.edu/tutorial/hubbard7.pdf
A very simple derivation for the current operator on a lattice is by interpreting the Heisenberg equation of motion as continuity equation.
If you are interested to the current of a quantity $\hat Q_i$ (a local operator defined for each site $i$, which could for example be particle number $\hat Q_i=\hat n_i=\hat a_i^\dagger \hat a_i$), you first derive the Heisenberg equation of motion $$ \frac{d}{dt}\hat Q_i=i[\hat H, \hat Q_i].$$ This already has the form of a continuity equation (i.e. $\dot\rho+\nabla\cdot\vec I=0$), so you can interpret the RHS as a finite difference version of $-\nabla\cdot\vec I$.
For example if you have $H=\sum_{<ij>}J_{ij}a_i^\dagger a_j+\text{H.c.}$, you find $$\frac{d}{dt}\hat n_i=i\sum_{j\in\{\text{nearest neigbours}\}}(J_{ij}a_i^\dagger a_j-J_{ij}^*a_j^\dagger a_i),$$ such that you can interpret the current from site $i$ to $j$ as $$\hat I_{ij}=i(J_{ij}a_i^\dagger a_j-J_{ij}^*a_j^\dagger a_i).$$
Totally revised
I had initially understood your question as asking about the continuum limit of the Hubbard model. Judging from your other question I now realize you meant to ask about current operators on a lattice model. As in, from the perspective of the lattice model alone (not thinking about any real space it might have been derived from) how do we talk about conserved current?
Note that since we are on a lattice there is no longer a normal notion of a derivative, so the continuity equation $\nabla\cdot j = \partial_t\rho$ no longer makes any sense. Also the Peierls/ minimal coupling prescription is ambiguous because it depends on the path we take from one site to another through the real continuous space. Before I go any further let me say it almost always sufficient to define a conserved current via an approximate formula like $$j(q) \approx e\int_{BZ} \!\!\!dk v(k)\psi^\dagger(k)\psi(k+q),$$ especially since the lattice model is itself usually an approximation, and we are almost always concerned about a very small portion of the Brillouin Zone.
But lets we would like a conserved current on the lattice. We would like to interpret the continuity equation as a finite difference, something like $j(x+a) - j(x) = \partial_t \rho$. I believe the correct way to talk about this is to imagine you current operator lives on the edges of your lattice. For definiteness let's say we have a (hyper)cubic lattice. Your vertices are labeled by a vector of integers $\vec{n}$, with the edges connecting nearest neighbors. Orient your edges so that they point in the direction of increasing coordinate. Now associated with each edge $e$ define an operator $j$. Our continuity equation is then $$\sum_{e\in in(n)}\!\!j(e) - \!\!\sum_{e\in out(n)}\!\!\!j(e) = \partial_t Q(\vec{n})$$
(which is just Kirchoff's Law). Now we want a local definition of $j$ in terms of the fields, which we should get because of something like Noether's theorem. For nearest neighbor hopping like in the Hubbard model we get $$j(e) = i\left[\bar{c}_{n_2}c_{n_1} - \bar{c}_{n_1} c_{n_2}\right]$$ where $\bar{c}$ and $c$ are the electron creation and annihilation operators respectively, and $n_1$ and $n_2$ are the two vertices connected by $e$. You can calculate by hand that this satisfies the continuity equation.
Note that we cannot in general just define $j$ to be a vector that lives on the vertices (like by averaging over the edges as below). This works in the the square lattice but not for example the triangular lattice in the plane. You can see this since there are six nearest neighbors but only two dimensions, so a vector does not have enough degrees of freedom to satisfy a continuity equation. Its also good form to have the current live on edges because this is the way to turn geometric objects on the continuum into geometric objects on a lattice: Scalars like the electric potential live on vertices, 1-forms like the current or E-field live on the edges, 2-forms like the magnetic field live on the plaquettes, etc... This definition also allows you to connect to the minimal substitution prescription, since the gauge fields lives on the edges.
To take the continuum limit of this formulation you simply define the vector $$\tilde{j}(x_n) = \sum_{e\in in(n)}\hat{e}j(e) - \sum_{e\in out(n)}\hat{e}j(e)$$, where $x_n$ is the real space position of the $\vec{n}$ site and $\hat{e}$ is the real space vector correpsonding to the lattice edge. If you look at this in fourier space you'll see that it recovers the usual continuity equation.
I mentioned Noether's theorem in the original version of this answer. There is, I've convinced myself, a statement that any local symmetry generator on a lattice has a conserved current, but its clunky and I don't really see any reason to state it in generality. Let me say this specific case: suppose your Hamiltonian can be written as $$H= \sum_n\mathcal{H}_1(\psi(n),n) + \sum_{<nm>}H_2(\psi^i_n,\psi^j_m;n,m)$$, where $\psi^i_n$ are fields that live on the $n$ vertex, and the sum $<mn>$ is taken over all pairs $<m,n>$ nearest neighbors. So you have a Hamiltonian which has term that involve at most nearest (no next-to nearest neighbors for example.) Now if you have a local symmetry generator $Q(n)$, which acts on the fields $\psi(n)$ and $[\sum_n Q(n), H]= 0$, then there is a conserved current
$$j_Q(e) = \frac{i}{2}\left[Q(n)-Q(m),H_2(\psi^i_n,\psi^j_m,n,m)\right]$$
where $n$ and $m$ are the two vertices $e$ connects. This $j_Q$ satisfies the continuity equation. You can just plug in and calculate. Note the $Q(n)$ has to be local; as in $[Q(n),\psi^i_n]$ depends only on the fields $\psi^j_n$ and $[Q(n),\psi^i_m] = 0$ when $n\neq m$. But this gives you, for example, the spin and charge currents for an arbitrary Hamiltonian as above.
You could also look at the continuity equation: If the particle density is conserved locally, then $\partial_t \rho(r) = -\nabla j(r)$, and also $\partial_t \rho(r) = -i [\rho(r),H]$. In the Hubbard model you will find that the interaction term commutes with the density operator, and therefore does not contribute to the current.
Another way would be to add the coupling to a vector potential (Peierls substitution, if you've heard of it), and take the derivative with respect to A: $j(r) \propto \frac{\delta H}{\delta A(r)}$
