I have to agree the exposition in this paper was quite sparse. I think answering these two questions requires a bit of review of previous work, but you can skip to the bottom if you just want the answers.
Recall the energy function of binary RBMs.
$$E(v,h) = v^TWh + v^Ta + h^Tb $$
In the binomial RBM work, they duplicate each visible unit a bunch of times, and correspondingly the bias $a_i$ of each visible unit and the rows $W_i$ corresponding to "connections" to the hidden units are also copied.
Now, let $W'$, $a'$, and $v'$ denote the appropriately duplicated values of $W$, $a$, and $v$. Critically, the energy function does not change, other than swapping out some symbols:
$$E(v',h) = v'^TW'h + v'^Ta' + h^Tb $$
From an implementation and computational standpoint, it's quite wasteful to do extra computation for all those units. You might ask if we can save any effort by representing $K$ duplicates of a binary unit using a single variable which just stores how many of those $K$ duplicates are set to 1 (and how many are set to 0). I'll use $v^*$ to denote the combined duplicates: in other words, $v^*$ is the sum of each group of $K$ duplicates in $v'$
In order to do so, we need to make sure we can still sample hidden units given visible units and vice versa, which is a prerequisite for running contrastive divergence.
Sampling hidden from visible is easy because $P(h=1|v') = \sigma(b+W'^Tv') = \sigma(b+Wv^*)$
Sampling visible from hidden is slightly tricker.
We have: $P(v'=1|h) = \sigma(a'+W'h)$. Now note that $v^*|h \sim \text{Binomial}(K, \sigma(a+Wh))$. By grouping out duplicated units together, we have gone from sampling from a bernoulli distribution to sampling from a binomial one.
Note that the gradient of the log probability:
$$\frac{\partial \log p(v)}{\partial W} = E_\text{data}[vh^T] - E_\text{model}[vh^T]$$
remains unchanged with our binomial units.
Now suppose after duplicating each unit $K$ times, we added different fixed biases to the duplicated units: the first duplicate bias is offset by $-0.5$, the second by $-1.5$, the third $-2.5$, and so on. In fact, why stop at $K$ duplicates -- why not make infinite duplicates of each unit. This is not as unreasonable as it might first sound, because each successive duplicate has a greater negative bias offset, so it has a vanishingly small chance of ever being switched on.
Whereas in the binomial case it was nice that we were able to group duplicate units together, here is it actually crucial, since computing with infinite binary units is not exactly possible. Whereas in the binomial case the sum of the duplicated units followed a binomial distribution, here it is possible to prove that the sum of the duplicated units approximately follows a distribution
$$\text{relu}(\mathcal{N}(x,\sigma(x)))$$
Take care to note that $\sigma$ here denotes the sigmoid function and not the standard deviation. So instead of sampling from a binomial distribution, we sample from this relu'd normal distribution, and everything else proceeds as previously described.
Now that we have an understanding of what exactly it means to have relu activations in an RBM, we can return to the questions at hand.
How should the energy function change?
We have seen above that it does not change at all!
How should max(0,x) be interpreted as a probability?
It shouldn't, and the relu RBM model never uses relu as a probability value of any sort. relu is only used to approximate the sampling of infinite binary units.
