Deduction of electric current density using field operators By obtaining the current of probability and multiplying by the electronic charge $e$, one can obtain the electric current density. By "replacing" the wavefunction by the field operator, one has
$$\mathbf{j}(\mathbf{r})=-i\frac{\hbar}{2m}\left[\Psi^{\dagger}(\mathbf{r})\nabla\Psi(\mathbf{r})-\left(\nabla\Psi^{\dagger}(\mathbf{r})\right)\Psi(\mathbf{r})\right].$$
How can this result be obtained by defining
$$\mathbf{j}(\mathbf{r}) = e\sum_i\mathbf{v}_i\delta(\mathbf{r}-\mathbf{r}_i)$$
Thanks
 A: The two operators are not identical, however they approach to each other in the semiclassical limit. To show that, we can compute matrix elements of the two expressions in a complete set of states and show that the matrix elements approach each other. (Strictly speaking this is weak convergence).
I'll perform the computation in one dimension, the generalization to higher dimensions is straightforward.
A possible basis is the coherent state basis. The coherent states are eigenstates of the field operator $\Psi(x)$.  The coherent states basis is overcomplete and especially convenient for semiclassical approximations. 
In this basis each Schrodinger particle is labeled by a position and momentum and  for $n$ Schroedinger particles, it has the form $| x_1, ..., x_n; p_1, ..., p_n \rangle $,  such that the action of the field operator on the coherent state:
$$\Psi(x) | x_1, ..., x_n; p_1, ..., p_n \rangle  = \prod_{i=1}^n \sqrt[\LARGE 4]{\frac {m \omega}{ \pi  \hbar }}  e^{ \large -\frac{m \omega}{2\hbar}(x-x_i)^2 + \frac{ixp_i}{\hbar}} | x_1, ..., x_n;  p_1, ..., p_n \rangle $$
The factors in the right hand side consist of the coherent state wave functions that Schrodinger wrote in 1926.
The constant parameter $\omega$ is a constant having units of angular frequency. 
The matrix element of the electric current ( I took the liberty of adding the particle electric charge to the expression) can be readily computed:
$$\langle  x_1, ..., x_n; p_1, ..., p_n |j(x) | x_1, ..., x_n; p_1, ..., p_n \rangle = e \sum_{i=1}^n \frac{p_i}{m} \sqrt{ \frac { m \omega}{  \pi  \hbar}} e^{-\frac{m \omega}{\hbar}(x-x_i)^2 }$$
The first factors are just the velocities:
$$v_i = \frac{p_i}{m}$$
The second term approaches a delta function in the semiclassical limit:
$$\lim_{\hbar \to 0} \sqrt{ \frac {m \omega }{  \pi  \hbar}} e^{-\frac{m \omega}{\hbar}(x-x_i)^2 } = \delta(x-x_i)$$
A: The answer given by @David Bar Moshe is enlightening but, seems to me, not so accurate since the $x_i$ and $p_i$ are merely labels of the coherent states, not the actual quantum operators (neither in first nor in second-quantized form).$\newcommand{\mb}{\mathbf}$
$\newcommand{\ima}{\mathrm{i}}$
$\newcommand{\di}{\mathrm{d}}$
Key point here is that, in first-quantization, $[\mb v, \delta(\mb x - \mb x')] \neq 0$, therefore $e\sum_i\mathbf{v}_i\delta(\mathbf{r}-\mathbf{r}_i)$ is not an Hermitian operator. You may hermitianize it by defining
$$\mb j(\mb x) = \frac{e}{2} \sum_i \left(\mb v_i \delta(\mb x - \mb x_i) + \delta(\mb x - \mb x_i) \mb v_i \right).$$
In coordinate representation, $\mb v $ is $-\frac{\ima \hbar}{m} \nabla$, and the second-quantized form can be obtained
$$
\begin{aligned}
        \hat j(\mb x) &= -\frac{e}{2m}\int \di \mb x' \, \left(\hat \psi^\dagger(\mb x') -\ima \nabla \left(\delta(\mb x - \mb x') \hat \psi(\mb x') \right) + \hat \psi^\dagger(\mb x') \delta(\mb x - \mb x')  -\ima \nabla \hat \psi(\mb x') \right) \\
        &= \frac{\ima e}{2m} \left(\hat \psi^\dagger(\mb x) \nabla \hat \psi(\mb x) - \nabla \hat \psi^\dagger(\mb x) \hat \psi(\mb x) \right)
    \end{aligned}
$$
(see for instance, Fetter, Alexander L., and John Dirk Walecka. Quantum theory of many-particle systems., if you do not know how these are rigorously derived)
