If $E$ (energy) and $P$ (momentum) only commute in constant potential, how could we have an $E$-$K$ diagram in a solid material?
$[E,p] \neq 0$. Then we cannot prepare electrons whose $E$ and $P$ are both very specified. However, in solid state material, we often have $E$-$K$ diagram, where every point is a pair of specified $E$ and $P$. The potential in a solid is definitely not a constant, so can somebody explain this?
In the band diagrams you refer to, $K$ is not actually momentum; instead, it is something called crystal momentum (also often 'quasimomentum'; for more details and its relationship to real momentum see this, this, this and this questions, among many others on this site).
In short, $\hat H$ does not commute with $\hat p$, but it does commute with a restricted subset of operators which are functions of $\hat p$. In particular, it can be seen that the operator $$ \hat U(a) = e^{i\hat p a/\hbar} $$ produces a translation by $a$ in position space by the displacement $a$. If the hamiltonian has the form $$ \hat H = \frac{1}{2m} \hat p^2 + V(\hat x) $$ where $V(x+a)=V(x)$ is a periodic potential, then generically $\hat H$ will not commute with arbitrary displacements, i.e. $[\hat H, e^{i\hat p \tilde a/\hbar}] \neq 0$ for arbitrary $\tilde a$, but it will commute with $e^{i\hat p \tilde a/\hbar}$ when $\tilde a = a$, i.e. when the displacement length matches the period of the potential.
Thus, we have $$ [\hat H, e^{i\hat p \tilde a/\hbar}] = 0, $$ which means that $\hat H$ and $e^{i\hat p \tilde a/\hbar}$ admit a common eigenbasis. Moreover, $\hat U(a) = e^{i\hat p a/\hbar}$ is a unitary operator, which means that all of its eigenvalues must have unit modulus and it is always possible to write them as $u = e^{ik a/\hbar}$, where $k$ is only ever defined up to a multiple of $2\pi\hbar/a$; to remove this ambiguity, we restrict ourselves only to $k$'s in the interval $(-\pi\hbar/a, \pi \hbar/a]$.
This $k$ is the crystal momentum that's used in the band diagrams you're bothered by. While it does have a nontrivial relationship with physical momentum, it is not the same thing.
This is not an answer to your question, but rather a response to the premise you wrote, to clear up a misconception: Energy and momentum can commute in a non-constant potential, it just depends what the potential is. You'd have to check it explicitly for a given potential by acting the momentum operator on the potential on a test function.
In general, in any closed system (so no external potential), momentum is conserved, which means we must have $$[H,p_i]=0$$
Why should they commute? A conserved quantity commutes with $H$ so that its time evolution is unaffected, thus it is conserved: $$P(t)=e^{-itH} P e^{itH}=P e^{-itH} e^{itH}=P$$
Edit: In the comments it was suggested that the claim above was false. So as an important example, below is a proof of this for the coulomb potential of two charge particles of mass $m$. Note that the conserved quantity is the summed (total) momentum of both particles, $\vec{p}_1+\vec{p}_2$. The individual momenta are not conserved and do not commute with $H$.
$$[H,p]$$ $$=[p_1^2/2m +p_2^2/2m+V(x_1,x_2), p_1 +p_2]$$
The first two terms of the hamiltonian commute with the total momentum so they vanish. $$=[V(x_1,x_2), p_1 + p_2]$$ $$=[V(x_1,x_2), p_1]+[V(x_1,x_2),p_2]$$
Consider the first commutator, for $p_1$. Use the notation $\partial_{x1}\equiv \frac{\partial}{\partial x_1}$.
$$[V(x_1,x_2), p_1]$$ $$=-i\hbar[V(x_1,x_2), \partial_{x1}]$$
where $V(x_1, x_2)=\frac{-k}{|x_2-x_1|}$
Let's find this commutator. Consider on some test function $\psi(x_1, x_2)$
$$\partial_{x1} (V(x_1, x_2) \psi)=V(x_1, x_2) \partial_{x1}\psi + \partial_{x1}V(x_1, x_2) \psi$$ $$\implies V(x_1, x_2) \partial_{x1}\psi-\partial_{x1} (V(x_1, x_2) \psi)=[V,\partial_{x1}]\psi=-\partial_{x1}V(x_1, x_2) \psi$$
The derivative of the potential is
$$\partial_{x1}V(x_1, x_2)=\frac{k(x_2-x_1)}{|x_2-x_1|^3}$$
Therefore $[V(x_1, x_2),\partial_{x1}]=\frac{k(x_2-x_1)}{|x_2-x_1|^3}$. This gives
$$[V(x_1, x_2),p_1]=-i\hbar\frac{k(x_2-x_1)}{|x_2-x_1|^3}$$ $$[V(x_1, x_2),p_2]=+i\hbar\frac{k(x_2-x_1)}{|x_2-x_1|^3}$$
..where the second commutator came from exchanging $x_1 \leftrightarrow x_2$. These commutators are negatives of each other, so adding them gives $$[V(x_1, x_2),p_1+p_2]=[V(x_1, x_2),p]=0$$ $$\implies [H,p]=0$$
This is a proof in one spacial dimension. The same proof works for the Coulomb potential in three spacial dimensions, where the generalized result is $[H,p_i]=0$ for all three components of the momentum.
