Evaluation of expectation values I will denote operators with hats. Suppose we got an operator of the form $i[\hat p, \tan^{-1}(e^{\hat x})]$ and we want to calculate the amplitude for a transition from a state $|p_i\rangle$ to the same state $|p_i\rangle$, as might happen in elastic scattering problems. Below are two evaluations which look like they are correct but take to two different results. Denoting $\hat f=\tan^{-1}(e^{\hat x})$,
1)
\begin{equation}
\begin{array}{lcl}
\langle p_i|i[\hat p,\hat f]|p_i\rangle&=&
i\langle p_i|(\hat p \hat f-\hat f\hat p)|p_i\rangle
=i\langle p_i|(p_i \hat f-\hat fp_i)|p_i\rangle\\
&=&i\,p_i\langle p_i|(\hat f-\hat f)|p_i\rangle=0
\end{array}
\end{equation} 
2)
\begin{equation}
\begin{array}{lcl}
\langle p_i|i[\hat p, \hat f]|p_i\rangle&=&i
\langle p_i|\left(-i\hbar\frac{\partial \hat f}{\partial x}\right)|p_i\rangle
=\hbar\langle p_i| \frac{e^{-\hat x}}{1+e^{-2\hat x}}|p_i\rangle\\
&=&\hbar\langle p_i| \int dx |x\rangle \langle x| \frac{e^{-\hat x}}{1+e^{-2\hat x}}|p_i\rangle
=\hbar\int dx  \frac{e^{- x}}{1+e^{-2 x}} \langle p_i| x\rangle \langle x |p_i\rangle\\
&=&\hbar\int dx \frac{e^{- x}}{1+e^{- 2x}} \frac{e^{\frac{i}{\hbar}p_i x}}{\sqrt{2\pi \hbar}}\frac{e^{-\frac{i}{\hbar}p_i x}}{\sqrt{2\pi \hbar}}
=\frac{1}{2\pi}\int_{-\infty}^{+\infty} dx \frac{e^{- x}}{1+e^{- 2x}}\\
&=&\frac{1}{2\pi}\,\frac{\pi}{2}=\frac{1}{4}
\end{array}
\end{equation} 
The integral in the evaluation 2) can be found in Jeffrey, Eq. (11), Sec. 15.3.1 .
I chose the form of $\hat f$ to bring the wanted results, but the question is general. Many other simpler forms can be employed to bring similar discordant results.
Moreover, the results may be generalized saying that: from the first evaluation the result is always zero, from the second evaluation the result is always $\propto f\big|_{x=+\infty}-f\big|_{x=-\infty}$.
What do I do wrong?
 A: The problem is ill posed from scratch because $|p_i\rangle$ is not an element of the Hilbert space and, a fortiori, it does not belong to the domain of $\hat{p}$. The same problem arises when considering $\hat{f}|p_i\rangle$.
As a matter of fact, rigorously speaking $\langle p_i|i[\hat p,\hat f]|p_i\rangle$ does not exist.
From the mathematical point of view the problem stops here since ex falso quodlibet.
However something can be said with a suitable interpretation of $\langle p_i|i[\hat p,\hat f]|p_i\rangle$. Naively but straightforwardly  one can give this interpretation, omitting inessential signs and constants,
$$\langle p_i|[\hat p,\hat f]|p_i\rangle= i\int_{\mathbb R} e^{-ip_ix} \left(\frac{d}{dx} f(x) e^{ip_ix}\right) dx - i\int_{\mathbb R} e^{-ip_ix} f(x)\frac{d}{dx} e^{ip_ix}  dx\:. \quad [1]$$
This integral can be computed and produces the result (2).
The result (1) cannot be obtained with this straightforward interpretation of the formalism, because it relies upon the formal self-adjointness  of $i\frac{d}{dx}$, that is the identity:
$$i \int_{\mathbb R} \psi(x) \frac{d}{dx} g(x) dx =  i\int_{\mathbb R} \left( \frac{d}{dx}\psi(x)\right)  g(x) dx \quad [2]$$
This identity, in fact holds for some classes of functions $\psi,g$ (in some cases also if $\psi, g$ do not belong to $L^2(\mathbb R)$).
Comparing with the first term in the right hand side of [1], we have that it should be,
$$g(x) = f(x) e^{ip_ix} \quad \psi(x) = e^{-ip_ix}\:.$$
However these functions have been fixed just to make false [2]!
A: It is not an answer directly to your case, but some relevant observations. Nevertheless it is too long to post it into comment.
We consider a similar situation, (I omit the hats)
$$\langle p | [x,p] | p \rangle  \tag{1} $$ 
There are two ways calculating expectation of (1). 
The first, $$ \langle p | xp -px | p \rangle = p ( \langle p|x|p \rangle - \langle p|x|p \rangle) =0 \tag{2} $$
The second, $$\langle p | i | p \rangle = i \delta(0) = i \infty \tag{3} $$
There is an inconsistency. The problem is, $$\langle p |x|p \rangle = \int \int dx dx' \langle p|x \rangle \langle x| x| x '\rangle \langle x'| p\rangle = \int x dx = \infty -\infty \tag{4} $$, which is ill-defined, unless we pick up principal value. 
We may make a more conservative calculation in the line of the first approach,
$$ \lim_{p'\rightarrow p} \langle p | [x,p] | p' \rangle =  \lim_{p'\rightarrow p} (p' - p) \langle p | x | p' \rangle = i  \lim_{p'\rightarrow p} (p' - p) \frac{ \partial}{\partial p} \delta(p-p') = i  \lim_{p'\rightarrow p} \delta (p - p') = i \delta(0) \tag{5} $$ 
Everything is consistent now.
I think the situation is similar for $f=ArcTan(e^x)$, though I did not work out the integral in $\lim p' \rightarrow p $ approach. 
A: (1) Looks wrong simply because you did not apply the operators onto the state $|p\rangle$ correctly. Operators act from right to left, so you should get:
$$
\left(\hat{p}\hat{f}-\hat{f}\hat{p}\right)|p\rangle=\hat{p}\hat{f}|p\rangle-\hat{f}\hat{p}|p\rangle 
 =\hat{p}\left(\hat{f}|p\rangle\right)-\hat{f}\left(\hat{p}|p\rangle\right) \tag{1}
$$
because $\hat{f}$ needs to act on $|p\rangle$ before you act $\hat{p}$ on it. The first term on the left of the second line (and neglecting constants) is really
$$
\hat{p}\left(\hat{f}|p\rangle\right)=\frac{\partial}{\partial x}\left( f|p\rangle\right)=\frac{\partial f}{\partial x}|p\rangle+fp|p\rangle\neq pf|p\rangle\tag{2}
$$
which is not what you have supposed.
The (2) looks to be correct because it is using the quantum analog of the Poisson bracket to define the commutator:
$$
[\hat{p},\,\hat{f}]=i\hbar\left(\frac{\partial\hat{p}}{\partial\hat{x}}\,\frac{\partial\hat{f}}{\partial\hat{p}}-\frac{\partial\hat{p}}{\partial\hat{p}}\frac{\partial\hat{f}}{\partial\hat{x}}\right)=-i\hbar\frac{\partial\hat{f}}{\partial x}\tag{3}
$$
because $\partial\hat{f}/\partial\hat{p}=0$.

EDIT 
The proper definition of the momentum operator is
$$
\hat{p}=-i\hbar\frac{\partial}{\partial x}
$$
Thus, by applying the missing constants to (2) and subtracting the $fp|p\rangle$ term from (1), we get
$$
-i\hbar\frac{\partial f}{\partial x}
$$
which is the same thing as equation (3). Thus, there is no paradox because you were not correctly applying your operators.

EDIT 2
The problem OP has is in the confusion about operator orders. When operators act on bras, we must take the (Hermitian) adjoint of the operators:
$$
\langle p|\hat p\hat f|p\rangle=\left(\langle p|\hat f^\dagger\hat p^\dagger\right)|p\rangle\neq \left(\langle p|\hat p\right)\hat f
$$
when using the middle term, a similar equation to my Equation (2) above will be obtained (i.e., ket replaced by a bra).
A: I think I found an answer that, at least for me, is satisfactory.
Transition amplitudes from a state $|a\rangle$ to the same state $|a\rangle$ make full sense if I am dealing with a bound state, whose wave-function is squared integrable. All problems I pointed out won't happen with bound (squared integrable) states. On the other hand, if I deal with a free state, whose quantum number is continuous, I'll have to insert a generally different state for the final state. I have pointed out scattering problems, but in scattering problems the scattered state $|p_ f\rangle$ is set to be in general different with respect to the initial state, although their energy is sometimes the same. 
One may be also tempted to consider $\langle p_i|V(\hat x)|p_i \rangle$ as expectation value of the operator $V(\hat x)$ onto the state $|p_i \rangle$. However, that looks wrong again if the state is not squared integrable, like $|p_i \rangle$. In fact, let us suppose we want to get the expectation value of the operator $\hat p$ onto the state $|p_i \rangle$. 
We naturally want the outcome to be $p_i$. This will not happen if we write the expectation value as 
$$
\langle p_i|\hat p|p_i \rangle=p_i\delta(0)=+\infty
$$
since the states are not squared integrable. The best we can get is by inserting a final different state and get the distribution 
$$
\mathcal{D}=\langle p|\hat p|p_i \rangle=p\delta(p-p_i)
$$
so that if integrated we get the wanted result:
$$
<p>=\int \mathcal{D} \,d\!p=p_i~,
$$
where $<p>$ means here "mean value for p".
In my opinion, when dealing with density matrices and traces in (x,p) quantum space, one must also do it with care, considering what above. 
In general, now I do not see anymore the need of considering elements of the type $\langle p_i|V(\hat x)|p_i \rangle$, whatever the potential is, if $|p_i\rangle$ is not squared integrable. In perturbation theory, which is what I mostly need, the first term of the Dyson series reads $\langle p_f|p_i \rangle$ and is set to zero, since we suppose $p_f \neq p_i$. If we don't suppose that, everything blows up to infinity at the zero order and that is what we do not want.
EDITED:
To emphasize what I mean, let me reconsider my equations with squared integrable states:
\begin{equation}
\langle p_i|x\rangle = \frac{e^{-\frac{i}{\hbar}p_i x}}{\sqrt{V}}\;,\;
\langle x|p_i\rangle = \frac{e^{\frac{i}{\hbar}p_i x}}{\sqrt{V}}\;,
\end{equation}
where $V$ is the integration 'Volume' in 1D (someone might want to call it $V^{1/3}$).
Equation (1) reads as it is, while Eq. (2) now reads:
\begin{equation}
\begin{array}{lcl}
\langle p_i|i[\hat p, \hat f]|p_i\rangle&=&
\hbar\int dx \frac{e^{- x}}{1+e^{- 2x}} \frac{e^{\frac{i}{\hbar}p_i x}}{\sqrt{V}}\frac{e^{-\frac{i}{\hbar}p_i x}}{\sqrt{V}}
=\frac{\hbar}{V}\,\int_{-\infty}^{+\infty} dx \frac{e^{- x}}{1+e^{- 2x}}\\
&=&\frac{\hbar}{V}\,\frac{\pi}{2}=0 ~,\qquad (1.1)
\end{array}
\end{equation} 
and the discrepancy has gone.
States in (1.1) are probably a good approximation of a Gaussian state with large spatial width, which is what we would get in usual experiments. I therefore expect that $\approx 0$ is the result we would find. As a matter of fact, a rigorous calculation with Gaussian states (which are squared integrable) would be interesting. We could send the spatial width to large values to see how the expectation value behaves. I'll see if I can do it in the following days.
