Tenga en cuenta que
$$\sum_{m=1}^\infty \frac{m}{q^{m+1}} = \frac{1}{(q-1)^2}, \qquad |q|>1.$$
Con $q= e^{nx}$, tenemos
$$f(x)=
\sum_{n=1}^\infty \sum_{m=1}^\infty \frac{m e^{-n (m+1) x}}{n}
=-\sum_{m=1}^\infty m \log\left(1-e^{-(1+m)x}\right). $$
A continuación, use el de Euler-Maclaurin fórmula:
$$f(x) \sim -\int_1^\infty \!dz\z \log\left(1-e^{-(1+z)x}\right)
- \overbrace{\frac{\log(1-e^{-2x})}{2}}^{\log x+ \mathcal{S}(1)}
+ \frac{1}{12} \overbrace{\left[ \frac{x}{e^{2x}-1}+\log \left(1-e^{-2 x}\right) \right]}^{\log x+\mathcal{S}(1)}
+ R$$
Con $$|R| \leq C \int_1^{\infty} \left| g''(z) \right|dz$$ where $g(z) = z \log\left(1-e^{-(1+z)x}\right) $. Numerical calculation shows $R = \mathcal{S}(1)$ aunque no tengo ninguna prueba todavía.
Así tenemos
$$f(x) \sim - \int_1^\infty \!dz\z \log\left(1-e^{-(1+z)x}\right)
= \frac{1}{x^2} \int_{2x}^\infty\!dx\,(x-y) \log(1-e^{-y}).$$
Necesitamos evaluar
$$\begin{align}\int_{2x}^\infty\!dy\,(x-y) \log(1-e^{-y})
&= x \underbrace{\int_{0}^\infty\!dy\, \log(1-e^{-y})}_{-\pi^2/6}
-\underbrace{\int_{0}^\infty\!dy\, y\log(1-e^{-y})}_{\zeta(3)}\\
&\quad+ \underbrace{\int_0^{2x} \!dy\,(y-x)\log(1-e^{-y})}_{\mathcal{O}(x^2)}\end{align} .$$
En conclusión, hemos
$$f(x) \sim \frac{\zeta(3)}{x^2}- \frac{\pi^2}{6x} - \frac{5}{12}\log x+ \mathcal{O}(1).$$