Even though Bayesian neural networks offer a promising framework for modeling uncertainty, active learning and incorporating prior physical knowledge, few applications of them can be found in the context of interatomic force modeling. One of the main challenges in their application to learning interatomic forces is the lack of suitable Monte Carlo Markov chain sampling algorithms for the posterior density, as the commonly used algorithms do not converge in a practical amount of time for many of the state-of-the-art architectures. As a response to this challenge, we introduce a new Monte Carlo Markov chain sampling algorithm in this paper which can circumvent the problems of the existing sampling methods. In addition, we introduce a new stochastic neural network model based on the NequIP architecture and demonstrate that, when combined with our novel sampling algorithm, we obtain predictions with state-of-the-art accuracy as well as a good measure of uncertainty.