A Riemannian stochastic representation of model uncertainties in molecular dynamics is proposed. The approach relies on a reduced-order model, the projection basis of which is randomized on a subset of the Stiefel manifold characterized by a set of linear constraints defining, e.g., Dirichlet boundary conditions in the physical space. We first show that these constraints are, indeed, preserved through Riemannian pushforward and pullback actions to, and from, the tangent space to the manifold at any admissible point. This fundamental property is subsequently exploited to derive a probabilistic model that leverages the multimodel nature of the atomistic setting. The proposed formulation offers several advantages, including a simple and interpretable low-dimensional parameterization, the ability to constraint the Fréchet mean on the manifold, and ease of implementation and propagation. The relevance of the proposed modeling framework is finally demonstrated on various applications including multiscale simulations on graphene-based systems.