Evaluating integral that include Fermi distribution function using residue theorem

Question

$$ I=\int_{-\infty}^{+\infty} \frac{d\epsilon}{2\pi} f(\epsilon) \left[ \frac{1}{(\epsilon-\epsilon_n +i0^+)^2(\epsilon-\epsilon_m +i0^+)} - \frac{1}{(\epsilon-\epsilon_n -i0^+)^2(\epsilon-\epsilon_m -i0^+)} \right] $$ According to the research article (Appendix D), if we use Cauchy's residue theorem, the integral becomes $$ I=- \left[ \frac{f(\epsilon_m)-f(\epsilon_n)}{(\epsilon_m-\epsilon_n)^2} + \frac{f'(\epsilon_n)}{(\epsilon_n-\epsilon_m)} \right] $$

here $f(\epsilon)=(e^{\beta\epsilon}+1)^{-1}$ and $0^+$ shows a infinitisimal small positive number. And $f'(\epsilon_n)=\frac{\partial f(\epsilon_n)}{\partial \epsilon_n} $.

My attempt:

Call the first term of $I$, $I_1$, and second term $I_2$. We have $$ I=I_1+I_2 $$ Calculating $I_1$:

To evaluate $I_1$ I look at the following contour integral: $$ I_u=\int_{C_u} \frac{dz}{2\pi} \left[ \frac{1}{(e^{\beta z}+1)(z-z_n )^2(z-z_m)} \right] $$ here $C_u=C_R+C$ is a semi-circle contour of radius $R$ in the upper plane of complex space, shown below

the two black dots are $z_{m(n)}=\epsilon_{m(n)}+i0^+$. In my understanding, this contour integral is equal to $I_1$ because the contribution of $C_R$ part of the contour $C_u$ goes to zero in the limit $R\to \infty $ (am I right?). And the remaining part contour $C$ becomes equal to $I_1$ in the limit $R\to \infty $. Next, I can use the residue theorem $$ \int_{C_u} g(z)dz=2\pi i \times \sum_k \text{Res}(g(z),z_k) $$ In our integral, we have three $z_k$ points,

1. $ z_1 = z_m $ is simple pole
2. $ z_2 = z_n $ is second order pole
3. $ z_3 = i\pi(2n+1)/\beta $ is a set of pole with $n \in \mathbb{Z} $

The residues at each $z_k$ are

1. $\text{Res}(g(z),z_1) = f(z_m)\frac{1}{(z_m-z_n )^2}= f(\epsilon_m)\frac{1}{(\epsilon_m-\epsilon_n )^2}$
2. $\text{Res}(g(z),z_2) = \lim_{z\to z_2} \frac{d}{dz}\left[f(z)\frac{1}{(z-z_m)}\right]= -f(\epsilon_n)\frac{1}{(\epsilon_n-\epsilon_m )^2} + f'(\epsilon_n)\frac{1}{(\epsilon_n-\epsilon_m )} $
3. $\text{Res}(g(z),z_3)$. I don't know how to calculate this residue. Any help?

All in all, we have $$ I_1 = I_u = i\left( \frac{f(\epsilon_m)-f(\epsilon_n)}{(\epsilon_m-\epsilon_n )^2} + \frac{f'(\epsilon_n)}{(\epsilon_n-\epsilon_m )} +\text{Res}\left(g(z),z_3\right)\right) $$

Calculating $I_2$:

Next, I do the same procedure for $I_2$ but with $z_{m(n)}=\epsilon_{m(n)}-i0^+$. The contour integral that is related to $I_2$ is $$ I_d=-\int_{C_d} \frac{dz}{2\pi} \left[ \frac{1}{(e^{\beta z}+1)(z-z_n )^2(z-z_m)} \right] $$ here $C_d$ is a semi-circle in the down plane of complex space. Now the poles are again at $z_m, z_n$, and $i\pi(2n+1)/\beta$. By residue theorem, I will get $$ I_2 = I_d = - i\left( \frac{f(\epsilon_m)-f(\epsilon_n)}{(\epsilon_m-\epsilon_n )^2} + \frac{f'(\epsilon_n)}{(\epsilon_n-\epsilon_m )} +\text{Res}\left(g(z),z_3\right)\right)=-I_1 $$

Hence, $I=0$

Can anyone please help me, why is my answer zero? It seems like I am double counting the poles.