The literature in MCMC and LBP never refer to the fact that the two methods look (on expectation) exactly the same. To illustrate, first consider a simple Ising model, that is, a graphical model where all the variables are boolean. We can denote these variables as $\{ X_i | X \in \{0,1\}, i \in (1,...,n)\}$ and we can write the joint distribution in some factored form as: $$ P(X) = \frac{1}{Z}\prod_{\alpha} \Psi_\alpha(X_\alpha) $$ Where $\alpha$ is some subset of the variables and $\Psi_\alpha$ is a potential function over these variables. Assuming we are interested in computing the marginal distribution of a given set of variables we could run Loopy BP to compute these or we could run a Gibbs sampler to simulate from $P$ to provide samples $\{X^{(1)},...,X^{(m)}\}$ and compute the marginal as: $X_i = \frac{1}{m}\sum_{j=1}^m X^{(j)}_i$.
For now, consider running an infinite number of Gibbs samplers at the same time, this would normally be computationally intractable since we would need to maintain an infinite number of sample sets. However, for discrete distributions we can represent these samples efficiently with a single vector. In our Ising model, for example, we can represent the sets of samples with a set of 2D vectors $\{V_i | V \in \mathbb{R}^2, i \in (1,...,n)\}$. The first component in the vector $V_i$ gives the proportion of the chains where $X_i = 1$ and the second component gives the proportion where $X_i = 0$ (of course we only actually need 1 dimension to represent this, but never mind). So, assuming we initialize the variables our Gibbs samplers uniformly then all of these vectors would equal $[0.5, 0.5]$. We can then continue with the regular Gibbs updates: loop through each of the variables (vectors) and re-sample them according to the neighbouring potentials (conditionals). After a burn in period these vectors will represent the set of samples which are drawn from the stationary distribution $P$, we can then compute the marginal distribution for each of the variables by averaging these vectors.
Now, other than the averaging step, this procedure is identical to loopy BP. And, in fact, the averaging is not likely to matter since these values tend to converge in mosts cases anyway.
If in fact, loopy BP is exactly the same as an "infinite" Gibbs sampler why is the literature on these two methods so different. All the Loopy BP analysis seems to be concerned with how it is minimizing the Bethe free energy while the Gibbs literature is focussed on mixing rates and the ergodicity of the Markov chain. Also, would simply averaging the messages in a Loopy BP inference procedure provide a correct estimate in cases where the updates are oscillating? Could the advancements in the Gibbs sampling literature, such as block Gibbs, be used in Loopy BP schemes to speed convergence? Lastly, if Loopy BP converges to an incorrect solution, does that imply that the Gibbs chain is non-ergodic?