Deriving the expectation of $[\hat X,\hat H]$ For a free particle of mass $m$, with Hamiltonian
$$\hat{H} = \frac {\hat{P}^2} {2m},$$
where
$$\hat{P} = -i \hbar \frac{\partial} {\partial x}.$$
The commutative relation is given by
$$[\hat{X}, \hat{H}] = \frac {i\hbar} {m} \hat{P}\tag{1}$$
In the common eigenstate of $\hat{H}$ and $\hat{P}$, $|e, p\rangle$,
can we do the following?
$$\langle e, p| [\hat{X}, \,\hat{H}] |e, p\rangle = \langle e, p|\hat{X} (\hat{H}|e, p\rangle) -  (\langle e, p|\hat{H}) \hat{X}|e, p\rangle 
\\ = \langle e, p|\hat{X} (e|e, p\rangle) -  (\langle e, p|e) \hat{X}|e, p\rangle \\
= e( \langle e, p|\hat{X}|e, p\rangle -  \langle e, p|\hat{X}|e, p\rangle ) \\ = 0 $$
Since the $\hat{H}$ is Hermitian, the above derivation doesn't seem to show any flaw.
Given the commutative relation, Eq (1), we know the result is wrong. What's wrong with the above derivation?
[EDIT]
Following the comment by Luboš Motl, I have worked out the solution and would like to share it here. The link provided by Qmechanic had the solution closely related to this question.
$$
\langle e', p'| [\hat{X}, \,\hat{H}] |e, p\rangle
\\ = \langle e', p'|\hat{X} (\hat{H}|e, p\rangle) -  (\langle e', p'|\hat{H}) \hat{X}|e, p\rangle 
\\ = (e - e') \langle e', p'|\hat{X}|e, p\rangle
$$
Note that: 
$$ e - e' = \frac{p^2}{2m} - \frac{p'^2}{2m} = \frac{(p+p')(p-p')}{2m} $$
$$ \langle e', p'|\hat{X}|e, p\rangle = -i\hbar \delta'(p - p') $$
where $\delta'(\cdot)$ is the derivative of the Dirac function, with respect to $p$.
Then we get
$$
(e - e') \langle e', p'|\hat{X}|e, p\rangle
\\ = -i\hbar \frac{(p+p')}{2m} \cdot (p - p')\delta'(p - p')
\\ = - \frac{i\hbar (p+p')}{2m} \cdot (-\delta(p - p'))
\\ = \frac{i\hbar (p+p')}{2m} \delta(p - p')
$$
As we take the limit $p \rightarrow p'$:
$$
lim_{p \rightarrow p'}
   \frac{i\hbar(p+p')}{2m} \delta(p - p')
\\ \rightarrow \frac{i\hbar}{m} p \delta(p - p')
$$
 A: As is customary in such question, I will point to this paper, which excellently discusses the problems the Dirac formalism has.
Now, in your concrete example, the problem lies in the energy/momentum states $| p_0 \rangle$ themselves, which are non-normalizable, since the wave function associated is the Fourier transform of $\delta(p-p_0)$, which means that $\psi_{|p_0\rangle}(x) = \mathrm{e}^{-\frac{ixp_0}{\hbar}}$. If you now try to calculate the inner product, you find: $$\langle p_0 | p_0 \rangle = \int_{-\infty}^{\infty}\psi_{|p_0\rangle}(x)\bar{\psi}_{|p_0\rangle}(x) \mathrm{d}x = \int_{-\infty}^{\infty} \mathrm{e}^{-\frac{ixp_0}{\hbar}} \mathrm{e}^{\frac{ixp_0}{\hbar}}\mathrm{d}x = \int_{-\infty}^\infty 1 \mathrm{d}x $$
Thus, momentum eigenstates are non-normalizable, and writing things like $\langle p_0 |X |p_0\rangle - \langle p_0 |X |p_0\rangle$ is really non-sensical, because you are subtracting two infinities. In particular, it is not $0$.
A: A tricky question, really. Apart from the fact that your $\lvert e,p\rangle$ vector does not belong to $L^2$ (hence you cannot take scalar products of it), I don't see any other flaw. That, in my opinion, means you have a nice argument to prove the following mathematical statement:

Let $\mathscr{H}$ be a separable Hilbert space, $0\neq z\in\mathbb{C}$. There are no self-adjoint operators $A$ and $B$ with non-empty discrete spectrum different from zero such that $[A,B]=z$.

Closely related to that fact, the following result of Von Neumann: up to multiplicity and unitary equivalence, the relations $[A,B]=i$ (in their exponentiated form) are uniquely realized by $A=x$ (multiplication operator) and $B=-i\nabla_x$, that indeed have no discrete spectrum.
EDITED (in reply to the comment, also the statement above has been edited slightly, to be more precise):
A number $\lambda\in \mathbb{R}$ is in the discrete spectrum of $A$ ( called $\sigma_{disc}(A)$ ) if there exists at least one $\psi_{\lambda}\in \mathscr{H}$ ( the Hilbert space, usually $L^2(\mathbb{R}^d)$ ) such that $$A\psi_\lambda=\lambda\psi_\lambda\; .$$ Suppose there exist $A$ and $B$ self-adjoint such that $0\neq \lambda \in \sigma_{disc}(B)$ and $[A,B]=z$ (on a suitable dense domain). Now it follows that (on another suitable domain)
$$[A,B^2]=2zB\; .$$
Let $\psi_\lambda\in \mathscr{H}$ be one of the eigenfunctions of $B$ associated to $\lambda$. On one hand,
$$2z\langle\psi_\lambda, B\, \psi_\lambda\rangle_{\mathscr{H}}=2z\lambda \lVert \psi_\lambda \rVert^2_{\mathscr{H}}\; ;$$
on the other
$$\langle\psi_\lambda, AB^2\, \psi_\lambda\rangle_{\mathscr{H}}- \langle\psi_\lambda, B^2A\, \psi_\lambda\rangle_{\mathscr{H}}=0$$
as you suggested. That is absurd, since $z$, $\lambda$ and $\lVert \psi_\lambda \rVert^2_{\mathscr{H}}$ are different from zero.
It follows that you cannot have two self-adjoint operators such that $[A,B]=z$ and $\sigma_{disc}(B)\neq \{0\},\emptyset$. The reasoning above does not work if there is no eigenfunction $\psi_\lambda\in \mathscr{H}$ (because with formal eigenfunctions you are not allowed to take scalar products or norms: they are not finite).
