Note that
$$\sum_{m=1}^\infty \frac{m}{q^{m+1}} = \frac{1}{(q-1)^2}, \qquad |q|>1.$$
With $q= e^{nx}$, we have
$$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). $$
Next, we use the Euler-Maclaurin formula:
$$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{O}(1)}
+ \frac{1}{12} \overbrace{\left[ \frac{x}{e^{2x}-1}+\log \left(1-e^{-2 x}\right) \right]}^{\log x+\mathcal{O}(1)}
+ R$$
With $$|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{O}(1)$ though I have no proof yet.
Thus we have
$$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}).$$
We need to evaluate
$$\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} .$$
In conclusion, we have
$$f(x) \sim \frac{\zeta(3)}{x^2}- \frac{\pi^2}{6x} - \frac{5}{12}\log x+ \mathcal{O}(1).$$