We consider the problem of time-stepping/sampling for molecular and meso-scale particle dynamics. The aim of the work is to derive numerical time-stepping methods that generate samples exactly from the desired target temperature distribution. The numerical methods proposed in this paper rely on the well-known splitting of stochastic thermostat equations into conservative and fluctuation-dissipation parts. We propose a methodology to derive numerical approximation to the fluctuation-dissipation part that exactly samples from the underlying Boltzmann distribution. Our methodology applies to Langevin dynamics as well as Dissipative Particle Dynamics and, more generally, to arbitrary position dependent fluctuation-dissipation terms. A Metropolis criterion is introduced to correct for numerical inconsistency in the conservative dynamics part of the model. Shadow energies are used to increase the acceptance rate under the Metropolis criterion. We call the newly proposed method meso-GSHMC.