Let $X_t$ be the (reflecting) diffusion process generated by $L:=Delta+ abla V$ on a complete connected Riemannian manifold $M$ possibly with a boundary $partial M$, where $Vin C^1(M)$ such that $mu(d x):= e^{V(x)}d x$ is a probability measure. We estimate the convergence rate for the empirical measure $mu_t:=frac 1 t int_0^t delta_{X_sd 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_1in 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 $dge 4$ and $ptoinfty$.