Not completely sure, you may want to check further: but I think $H^M$ is not directly the real Hamiltonian --- in your reference "construct the operator" means it's a newly-defined operator with $|\Psi_k\rangle$ and $\rho_k$; and $\beta$ is also newly-defined there, not necessary meaning the temperature.
I guess the reason you take it as the Hamiltonian is that the final "partition function": $Z=\text{Tr}\big[e^{-\beta H^M}\big]$ looks really like the partition function for canonical ensemble. Then the question now is why these two things are "the same"? To answer it, let us forget about definition $H^M$ for a second, which is directly related to $\rho$'s eigensystems $|\psi_k\rangle$ and $E_k^M$, and go back to the definition in quantum statstics of canonical ensemble.
Recall that, by definition (actually equal-distribution-theorem), when equilibrium is reached, the probability of the system being in an energy window around $E_{k}$ is proportional to $e^{-\beta E_k}$, but be careful here $E_k$ is the Hamiltonian $H$'s eigenvalue. Now we ask the question: what is the density matrix of such an equilibrium ensemble? The result turns out to be:
\begin{align}
\rho_{eq} \propto \sum_k e^{-\beta E_k}| \phi_k\rangle\langle \phi_k|
\end{align}
"$_{eq}$" means this is for equilibrium . Remember, again, here we only use real Hamiltonian $H$, and $|\phi_k\rangle$ and $E_k$ are corresponding eigensystem of it. Therefore, our logic is that: from original Hamiltonian, we construct the equilibrium density matrix $\rho_{eq}$ which satisfies the probability distribution desired. In this way, for a canonical ensemble, the equilibrium density matrix is directly diagonalized in Hamiltonian's eigen-basis, by definition. Then, like you said, that means:
\begin{align}
i\frac{d}{dt}\rho_{eq} = \big[H, \rho_{eq}\big] = 0
\end{align}
and $\rho_{eq}$ doesn't evolve any more -- which makes sense since it's equilibrium.
Now, what about the logic in your reference? It turns out it's constructed in an inverse way: from an arbitrary density matrix $\rho$, construct an operator $H^M$, which is not necessarily the real Hamiltonian of the system, and also define a $\beta$ from $\rho_k$, which is not necessarily the temperature. So what's their meaning? Well, if you want, you could say if you design a system with Hamiltonian $H^M$ and temperature $\beta$, then the given density matrix $\rho$ happens to be the equilibrium density matrix $\rho_{eq}$:
\begin{align}
\rho_{eq}\big(H^M,\beta\big) = \rho
\end{align}
In a short word, this is not surprising: they just design a system with equilibrium being the given $\rho$. And it's not contradict to the von Neumann equation, since in von Neumann equation $\rho$ and $H$ are real density matrix and real Hamiltonian of the system, other than designed from each other --- which is the logic of canonical ensembles or the logic in your reference.
