I think I solved it; unfortunately I don't have time to do all the manipulations after the integral is computed.
Let's change to spherical coordinates
$$I=-\pi(\vec{p}_m\cdot\vec{p}_n)\int \text{d}q\text{d}\phi \sin(\phi) \text{d}\theta\frac{1}{q}\frac{1}{(E_n-p_n \cos(\phi))(E_m-p_m \cos(\phi))}$$
It follows that:
$$I=-\pi(\vec{p}_m\cdot\vec{p}_n)\ln\left(\frac{\Lambda}{\lambda}\right)2\pi \int_0^{\pi}\text{d}\phi \frac{\sin(\phi)}{(E_n-p_n \cos(\phi))(E_m-p_m \cos(\phi))}$$
Now, let's change variables
$$x=\cos(\phi) \\ \text{d}x=-\sin(\phi)\text{d}\phi$$
so we have
$$I=-2\pi^2(\vec{p}_m\cdot\vec{p}_n)\ln\left(\frac{\Lambda}{\lambda}\right) \int_{-1}^{1}\text{d}x\frac{1}{(E_n-p_n x)(E_m-p_m x)}$$
Now, we have to make the following step:
$$\frac{A}{(E_n-p_n x)}+\frac{B}{(E_m-p_m x)}$$
you find that
$$A=\frac{-p_n}{E_n p_m-E_m p_n} \\ B=\frac{p_m}{E_np_m-E_mp_n}$$
Then, it's easy to see that
$$I=2\pi^2(\vec{p}_m\cdot\vec{p}_n)\ln\left(\frac{\Lambda}{\lambda}\right)\left[\frac{A}{p_n}\ln\left(E_n-p_nx\right)+\frac{B}{p_m}\ln\left(E_m-p_mx\right)\right]^1_{-1}$$
by making $A$ and $B$ explicit we get:
$$ I=2\pi^2\frac{(\vec{p}_m\cdot\vec{p}_n)}{E_np_m-E_np_n}\ln\left(\frac{\Lambda}{\lambda}\right) \left[-\ln\left(E_n-p_nx\right)+\ln\left(E_m-p_mx\right)\right]^1_{-1}$$
therefore
$$I=2\pi^2\frac{(\vec{p}_m\cdot\vec{p}_n)}{E_np_m-E_np_n}\ln\left(\frac{\Lambda}{\lambda}\right)\left[\ln\left(\frac{E_m-p_m}{E_n-p_n}\right)+\ln\left(\frac{E_n+p_n}{E_m+p_m}\right)\right]$$
and then it's should just be rearranging the factors in a convenient way. I hope that helped!