Potential of mean force (PMF) calculations along a reaction coordinate require exhaustive sampling which imposes severe computational limitations for any algorithm applied. The Expanded Ensemble Density of States method (EXEDOS) constitutes an elegant formalism to sample continuously and uniformly the reaction coordinate and extract the desired PMF. While the algorithm is trivial to parallelize, the iterative estimation of the required weights can be rendered exceedingly slow through insufficient sampling of the remaining degrees of freedom. In the present study, we propose a combination of Monte Carlo methods that can overcome the computational barriers present in the algorithm in a variety of systems where the reaction coordinate is the separation distance between two molecules. We applied the EXEDOS formalism in a simple well-studied system of two large colloids immersed in a solution of small particles. We demonstrate that calculations at all densities can be significantly accelerated through selective sampling of translational degrees of freedom in the vicinity of the large colloids. At low densities, application of rejection free Monte Carlo moves for the large particles can result in further improvement by almost an order of magnitude in computational time. At high densities where simple Monte Carlo simulations are prohibitively demanding and rejection-free moves fail, we explore recently proposed algorithms for reversible mapping between local minima that could provide a unique route for accurate PMF calculations in dense molecular systems.