Let $X_t$ be the (reflecting) diffusion process generated by $L:=Δ+\nabla V$ on a complete connected Riemannian manifold $M$ possibly with a boundary $\partial M$, where $V\in C^1(M)$ such that $μ(d x):= e^{V(x)}d x$ is a probability measure. We estimate the convergence rate for the empirical measure $μ_t:=\frac 1 t \int_0^t δ_{X_s\d s$ under the Wasserstein distance. As a typical example, when $M=\mathbb R^d$ and $V(x)= c_1- c_2 |x|^p$ for some constants $c_1\in \mathbb R, c_2>0$ and $p>1$, the explicit upper and lower bounds are present for the convergence rate, which are of sharp order when either $d<\frac{4(p-1)}p$ or $d\ge 4$ and $p\to\infty$.