I can see it is right by using some knowledge of the Gamma function. We have
$$ \Gamma(\alpha ) = \int_0^\infty e^{-x} x^{\alpha - 1 } dx . $$
Differentiating with respect to $\alpha$, we get
$$ \frac{d\Gamma}{d \alpha} = \int_0^\infty e^{-x} x^{\alpha - 1 }\ln x dx . $$
Setting $\alpha = 1 $, we get
$$ \int_0^\infty e^{-x} \ln x dx = \frac{d\Gamma}{d\alpha }\bigg|_{\alpha = 1 }= - \gamma. $$
But the knowledge $d \Gamma/d \alpha |_{\alpha = 1}=-\gamma$ is a mystery to me.
Can anyone find an elementary proof?
You can prove that $$\lim I_n=\int_0^\infty e^{-x}\ln x\,dx$$ using the dominated convergence theorem, where $$I_n=\int_0^n\left(1-\frac xn\right)^n\ln x\,dx.$$ Now substitute $y=1-x/n$. Then \begin{align} I_n&=n\int_0^1y^n\ln(n(1-y))\,dy\\ &=n\ln n\int_0^1 y^n\,dy+n\int_0^1y^n\ln(1-y)\,dy\\ &=\frac n{n+1}\left(\ln n-\int_0^1\frac{1-y^{n+1}}{1-y}\,dy\right)\\ &=\frac n{n+1}\left(\ln n-\sum_{k=1}^{n+1}\frac1k\right) \end{align} where we have integrated by parts along the way. This tends to $-\gamma$.