I guess for the modified Bessel funcion $K_0(z)$, $$\sum_{n=1}^\infty K_0(s\, n) \sim \frac{-2\log 2 - \log \pi + \gamma}{2} + \frac{\log s}{2} + \frac{\pi}{2\, s}, \quad s\to 0^+,$$ if taking $$\sum_{n=1}^\infty 1 =\zeta(0),$$ and $$\sum_{n=1}^\infty \log n = \log \prod_{n=1}^\infty n=\log \sqrt{2\pi},$$ I can get all the rest terms, but where does the term $\frac{\pi}{2\,s}$ come from?

The term $\frac{\pi}{2s}$ comes from the pole of $\zeta(s)$. Let's use the fact that Mellin transform of $K_0(s)$ equals $$ \int_0^{+\infty} K_0(s)s^{t-1}ds=2^{t-2}\Gamma^2(t/2). $$

From this we get that if $f(s)=\sum_{n} K_0(sn)$ then for $\mathrm{Re}\,s>1$ we have

$$ \int_0^{+\infty} f(s)s^{t-1}ds=2^{t-2}\Gamma^2(t/2)\zeta(t) $$

Therefore, by Mellin inversion we get

$$ f(s)=\frac{1}{2\pi i}\int_{3/2-i\infty}^{3/2+i\infty} 2^{t-2}\Gamma^2(t/2)\zeta(t)s^{-t}dt. $$

Now, if you move your contour far to the left, you will encounter poles of the integrand at the points $t=1,0,-2,-4,-6\ldots$. Using the fact that $\Gamma$ decays exponentially on the vertical lines and $\zeta$ grows at most polynomially, we obtain for any $N$

$$ f(s)=\mathrm{Res}_{t=1} 2^{t-2}\zeta(t)\Gamma^2(t/2)s^{-t}+\sum_{n=0}^{N} \mathrm{Res}_{t=-2n} 2^{t-2}\zeta(t)\Gamma^2(t/2)s^{-t}+O(s^{-1-2N}) $$

Computation of the first term reveals the desired $\frac{\pi}{2s}$ term, as $\Gamma(1/2)^2=\pi$.