A new generalized hyper-parallel tempering Monte Carlo simulation method is presented. The method is particularly useful for simulation of many-molecule complex systems, where rough energy landscapes and inherently long characteristic relaxation times can pose formidable obstacles to effective sampling of relevant regions of configuration space. In this paper, we demonstrate the effectiveness of the new method by implementing it in a grand canonical ensemble for the Lennard-Jones fluid and the restricted primitive model. Coexistence curves and critical behavior have been explored by the new method. Our numerical results indicate that the new algorithm can be orders of magnitude more efficient than previously available techniques.