NDSolve::ndnum: Encountered non-numerical value for a derivative at t==0.

σ = 1
rc = -4
λ = 1
sink[r_, λ_] = (3 λ)/(4 π r^6)
U[r_] = rc/r
g[r_] = Exp[-U[r]]
Dif = 1
h[r_] = 1
r0 = 1;
tmax = 0.01;
rmax = 1000 Sqrt[6 Dif tmax] + r0;

sol = NDSolve[{D[ρ[r, t], {t, 1}] ==
Dif r^-2 D[
   r^2 h[r] Exp[-U[r]] D[Exp[U[r]] ρ[r, t], {r, 1}], {r,
    1}] - (λ/g[σ] DiracDelta[r - σ]/(
     4 π σ^2) + sink[r, λ]) ρ[r,
   t], ρ[r, 0] == Exp[-(U[r] - rc/rmax)], ρ[rmax, t] ==
1, 4 π Dif h[σ] (
  D[Exp[U[r]] ρ[r, t], {r, 1}] /.
   r -> σ) == λ Exp[U[σ]] ρ[σ,
  t]}, ρ, {r, σ, rmax}, {t, 0, tmax},
  PrecisionGoal -> 4, StartingStepSize -> 0.0001,
  Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
  "MinPoints" -> 100, "MaxPoints" -> 200}}]

ρn = First[ρ /. sol]
Plot3D[ρn[r, t], {r, σ, rmax/100}, {t, 0, tmax}, PlotRange -> All]

appro = With[{k = 1}, Erf[k #]/2 + 1/2 &];

Plot[appro'[x - σ], {x, -σ, rmax}, PlotRange -> All, Axes -> None]

σ = 1;
rc = -4;
λ = 1;
sink[r_, λ_] = (3 λ)/(4 π r^6);
U[r_] = rc/r;
g[r_] = Exp[-U[r]];
Dif = 1;
h[r_] = 1;
r0 = 1;
tmax = 0.01;
rmax = 1000 Sqrt[6 Dif tmax] + r0;

mol[n_Integer, o_:"Pseudospectral"] := {"MethodOfLines",
  "SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> n,
    "MinPoints" -> n, "DifferenceOrder" -> o}}

mol[tf:False|True,sf_:Automatic]:={"MethodOfLines",
"DifferentiateBoundaryConditions"->{tf,"ScaleFactor"->sf}}

sol = With[{coef = 10^6},
   NDSolve[{D[ρ[r, t], {t, 1}] ==
       Dif r^-2 D[
          r^2 h[r] Exp[-U[r]] D[Exp[U[r]] ρ[r, t], {r, 1}], {r,
           1}] - (λ/
             g[σ] coef  DiracDelta[r - σ]/(4 π σ^2) +
           sink[r, λ]) ρ[r, t], ρ[r, 0] ==
       Exp[-(U[r] - rc/rmax)], ρ[rmax, t] == 1,
      4 π Dif h[σ] (D[Exp[U[r]] ρ[r, t], {r, 1}] /.
          r -> σ) == λ Exp[U[σ]] ρ[σ, t]} /.
     DiracDelta -> appro', ρ, {r, σ, rmax}, {t, 0, tmax},
    Method -> Union[mol[True, 100], mol[4000, 4]]]];

ρn = First[ρ /. sol];
Plot3D[ρn[r, t], {r, σ, rmax/100}, {t, 0, tmax}, PlotRange -> All]