Our proposed framework allows for the utilization of any number of Gaussian functions in a linear combination to approximate complex potentials. Its ability to represent various phases, including coexistence, while maintaining the principles of HhPF, makes it a promising approach within mesoscale molecular simulation techniques. This advancement opens up new possibilities for studying complex systems and their phase behaviors at the mesoscopic level.