The Problem
I was trying to show that the momentum operator is given by $\sum_r c^{\dagger}_r (-i\hbar)\nabla c_r$.This is how I tried showing it
$$\hat{P}=\sum_{r,r'}c_r^{\dagger} \langle r|\hat{P}|r'\rangle c_{r'}=\sum_{r,r'}c_r^{\dagger}(i\hbar\frac{d\delta(r-r')}{dr})c_{r'}=\sum_{r,r'}-i\hbar \frac{dc_r^{\dagger}}{dr}c_{r'}\delta(r-r') = \sum_r \nabla(c_r^{\dagger})(-i\hbar)c_r$$
As you can see, something has gone wrong. The formula derived doesn't even look hermitian! ($(\nabla(c_r^{\dagger})c_r)^{\dagger}=c_r^{\dagger}c_r\nabla)$. What went wrong here?
My Thoughts
I believe the issue here is with the behaviour of the matrix elements $\langle r|\hat{P}|r'\rangle =i\hbar\frac{d\delta(r-r')}{dr}$. Is there some shenanigans happening in there due to the derivative of the delta function? If the matrix element had instead $\frac{d}{dr'}\delta(r-r')$, then I believe things would have been fine. However, did I make a mistake in evaluating this?